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■^ I Abstract. We study the Bose-Hubbard and Fermi-Hubbard model in the limit of 

large coordination numbers Z (i.e., many tunnelling partners). Via a eontrolled 
expansion into powers of l/Z , we establish a hierarchy of correlations, which facilitates 
(-H ' an approximate analytic solution of the quantum evolution. For the Bose-Hubbard 

0^. model, we derive the growth of phase coherence after a quench from the Mott to the 

superfluid phase. For a quench within the Mott phase, we find that various local 
observables approach a quasi-equilibrium state after a finite period of time. However, 
r^ , this state is not thermal, i.e., real thermalisation - if it occurs - requires much longer 

O^' time scales. For a tilted lattice in the Mott state, we calculate the tunnelling probability 

and find a remarkable analogy to the Sauter-Schwinger effect (i.e., electron-positron 
CS| , pair creation out of the vacuum due to a strong electric field). These analytical results 

^ ' are compared to numerical simulations for finite lattices in one and two dimensions 

and we find qualitative agreement. Finally, we generalise these studies to the more 
involved case of the Fermi-Hubbard model. 
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1. Introduction 

The theoretical description of strongly interacting many-particle quantum systems in 
condensed matter physics is a difficult undertaking in general. Even in those cases where 
it is possible to derive a Hamiltonian which contains all relevant parts for a complete 
description of the many-particle system under consideration, a complete solution of the 
dynamics or an exact evaluation of the ground state is often out of reach. 

However, essential features such as phase transitions or quasi-particle excitations are 
already contained in drastically simplified models. One of the most studied models which 
describes strongly interacting electrons in a solid is the Fermi-Hubbard model [H El [3] . 
Although it describes only the physics of a single energy band without long-ranged 
Coulomb-interactions, it is believed to exhibit various interesting phenomena such as 
the metal-insulator transition, anti-ferromagnetism, ferrimagnetism, ferromagnetism, 
Tomonaga-Luttinger liquid, and superconductivity [1]. The Fermi-Hubbard model was 
ffist established by J. Hubbard in order to describe correlations of electrons in narrow 
energy bands [H HI [3]. It connects the Heitler-London theory of strongly interacting 
electrons on one side and the band theory, valid for weak interactions, on the other side. 

Apart from electrons in solids, there are also other physical systems which can be 
described (within suitable approximations) by the Fermi-Hubbard model, for example 
ultra-cold atoms in optical lattices. Since optical lattice systems allow for a tuning 
of the model parameters over a wide range (in contrast to most condensed matter 
systems), they are predestinated for the study of the fermionic Hubbard model. This 
has been experimentally realised by trapping fermionic Potassium atoms in an optical 
lattice O [HI [7] . By increasing the on-site interaction among the atoms, the transition 
from a metallic phase to a Mott insulating phase was deduced from compressibility 
measurements and in situ imaging [7]. 

After replacing the fermionic by bosonic atoms, the optical lattice system 
corresponds to the Bose- Hubbard model, which describes interacting bosons in periodic 
potentials. The Bose-Hubbard model was ffist motivated by experimental realisations 
such as Helium-4 absorbed in porous media or Cooper pairs in granular media [8]. 
Within the last decade, the Bose-Hubbard model has been attracting increasing 
attention due to its experimental realisation with interacting bosons, for example 
Rubidium atoms, in optical lattices O [10]. The phase diagram, which is much 
better understood than in the fermionic case, contains (for vanishing disorder) a 
superfiuid regime and Mott insulating phases. The transition between the two phases is 
characterised by a natural order parameter, such as the superfiuid density. This second- 
order quantum phase transition, which results from the competition between kinetic 
energy and on-site interaction, has been observed for bosons in optical lattices [9]. 

Though at ffist glance seemingly very simple, the fermionic as well as the bosonic 
lattice models are not exactly solvable in general, the only exceptions being the Fermi- 
Hubbard model in one dimension [HI [12] and trivial situations, such as the case of 
vanishing or infinitely strong interaction, or an empty lattice, for example. Thus, 
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various analytical and numerical techniques have been developed to study their physical 
properties. Widely used numerical methods, ignoring distance-dependent quantum 
correlations between different lattice sites, are mean-field approaches such as the 
dynamical mean-field theory (DMFT) [131 [H]- Exact diagonalisation (analytical or 
numerical) is only possible for small system sizes [T5l [16]. Monte Carlo methods can 
improve the situation [17], but they can sample only a small part of the whole Hilbert 
space and have problems with general time-dependent situations. 

An important and widely used analytical approach is the Gutzwiller ansatz [18] - 
for the bosonic and the fermionic case. However, a serious drawback of this Gutzwiller 
mean-field approach is again the neglect of distance-dependent quantum correlations 
between different lattice sites. This simplification leads sometimes to unphysical 
behaviour - for example, the Mott-insulator state would not react to an external force 
within the Gutzwiller mean-field approach. 

The main goal of the present work is to study these non-local quantum correlations 
in the Bose and Fermi-Hubbard models. To this end, we develop an analytical expansion 
in powers of the inverse coordination number Z (i.e., the number of tunnelling partners 
of a given lattice site). For comparison, we also calculated these non-local quantum 
correlations by means of exact (numerical) diagonalisation for one and two-dimensional 
lattices. Although exact diagonalisation is possible only for small systems, the results 
are in good agreement with those obtained by the analytical expansion and the Density 
Matrix Renormalisation Group technique (DMRG) [T9| [20]. 

The paper is organised as follows. After a brief introduction into the Bose- 
Hubbard model in Sec. [21 we present our analytical method for general Hubbard type 
Hamiltonians in Section [3l With this method, the ground state properties of the bosonic 
Mott insulator phase are derived in Sec. [H Taking this Mott insulator as the initial state, 
the temporal evolution of the correlations after a quantum quench to the superfiuid 
regime as well as within the Mott phase are studied in Sections [5] and [6l In Section [3, we 
investigate particle-hole pair creation out of the (bosonic) Mott state induced by a weak 
tilt of the lattice and establish a quantitative analogy to electron-positron pair creation 
due to a strong electric field (Sauter-Schwinger effect) in Sec. [81 The results of Sees. [IHH 
are derived in first order 1/Z (where Z is the coordination number). In Section [HI we 
show how to extend our analytical method to higher orders in 1/Z. Subsequently, we 
compare our analytical results with exact numerical studies of finite bosonic lattices in 
Section [To] and find qualitative agreement. 

In the second part of the paper, we consider similar problems for fermions. After 
a brief introduction into the Fermi-Hubbard model in Section [Til we discuss its ground 
state correlations (in the fermionic Mott state) and quench dynamics Sees. [12] and [131 
Section [m is then devoted to particle-hole pair creation induced by a weak lattice tilt. 
Finally, we address resonant tunnelling in the Bose and Fermi-Hubbard model due to a 
large lattice tilt in Section [TS] 
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2. Bose-Hubbard Model 

The Bose-Hubbard model is one of the most simple and yet non-trivial models in 
condensed matter theory [8], [21], [22] . It describes identical bosons hopping on a lattice 
with the tunnelling rate J. In addition, two (or more) bosons at the same lattice site 
repel each other with the interaction energy U. The Hamiltonian reads 

H = -^Y.^,ubik + ^Y.^,in,-l). (1) 

Here 6j^ and bi, are the creation and annihilation operators at the lattice sites /i and u, 
respectively, which obey the usual commutation relations 



bu,bi 



^l-iiy 5 



blK 



bu,b. 



0. (2) 



The lattice structure is encoded in the adjacency matrix T^^^ which equals unity if fi and 
z/ are tunnelling neighbours (i.e., if a particle can hop from /i to u) and zero otherwise. 
The number of tunnelling neighbours at a given site /i yields the coordination number 
Z = J2u'^f^i' (^^ assume a translationally invariant lattice). Finally, n^ = b^^b^ is the 
number operator and we assume unit filling (n^) = 1 in the following. Note that the 
total particle number A^ = ^„ ^^ is conserved [H, N] = 0. 

The Bose-Hubbard model is considered [23] one of the prototypical examples for a 
quantum phase transition. If the interaction term dominates U ^ J, the bosons are 
pinned to their lattice sites and we have the Mott insulator state 

I^Mott),=o = (8)|l)M = n^i|0) ^ ^|^Mott)j=o = 0, (3) 

which is fully localised. If the hopping rate dominates U <^ J, on the other hand, the 
particles can propagate freely across the lattice and become completely delocalised 

i'p..p.„„,.)„.„ = -^ (i>u) 10) = -^ E *;. 10) . w 



which is the superfluid phase. Obviously, the Mott state (|3]) does not have any 
correlations, for example {bjj^b^)Mou = ^fiu, whereas the superfluid state in (jlj) shows 
correlations across the whole lattice (&t&;/)supcrfiuid = 1- Furthermore, the Mott insulator 
state is separated by a finite energy gap from the lowest excited state, while the 
superfluid state possesses sound-like modes with arbitrarily low energies (for an infinitely 
large lattice A^ f oo). Finally, the Bose-Hubbard model can be realised experimentally 
(to a very good approximation) with ultra-cold atoms in optical lattices [23], [251 [26] and 
it was even possible to observe the aforementioned phase transition in these systems ^. 
In spite of its simplicity, the Bose-Hubbard model ([1]) cannot be solved analytically. 
Numerical simulations are limited to reduced sub-spaces or small systems sizes, see 
Section [TO] below. Analytical approaches are based on suitable approximations. In order 
to control the error of these approximations, they should be based on an expansion in 
term of some large or small control parameter. For the Bose-Hubbard model ([I]), one 
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could consider the limit of large (n^) ^ 1 or small (n^) ^ 1 filling [271 I2B], for example, 
or the limit of weak coupling f/ ^ J or strong coupling U ^ J [291 EQl [31] . However, 
none of these limits is particularly well suited for studying the Mott-superfiuid phase 
transition. To this end, we consider the limit Z ^ 1 in the following and employ an 
expansion into powers of 1/Z as small control parameter. Note that an expansion in 
powers of 1/Z was also used to derive bosonic dynamical mean- field equations (which 
were then solved numerically) in [321 |33l [31] . 

3. Hierarchy of Correlations 

Let us consider general Hamiltonians of the form 

^ = |E^mv + E^m, (5) 

which includes the Bose-Hubbard model ([ID as a special case. The quantum evolution 
of the density operator p describing the state of the full lattice can be written as 

idtp = H,p =^Y1 [^^'^' p]^Yl [^A*' P 

fJ,U fl 

= -^^^u'^P + ^^^^P^ (6) 

fJ,U fl 

where we have introduced the Liouville super-operators £^jy and £^ as short-hand 
notation. As the next step, we introduce the reduced density matrices for one or more 
lattice sites via averaging (tracing) over all other sites 

Pf^u = tr^^{p} , (7) 

and so on. Note that tr{p} = 1 implies tr^{p^} = 1 and ti^uiPfiu} = 1 etc. Since we are 
interested in the (quantum) correlations, we separate the correlated and uncorrelated 
parts of the reduced density matrices via 

PnuX = P'-^TX + P™"PA + p';rPu + PTP,^ + WuPX , (8) 

and analogously for more lattice sites. As a consequence, we obtain from ([6]) the 
evolution equation for the on-site density matrix 

'dtp, = ^Y.^^A ^Up7^' + P^P^) } + ^mPm , (9) 

where £j^j, = C^y+ C^fi denotes the symmetrised form. Obviously, solving this equation 
exactly requires knowledge of the two-point correlation p^°". The time-evolution of this 
quantity can also obtained from (l6|) and reads 

tdtp^;:^ = C.p^Z^ + 1 £,.(p-" + p,p^) - f tr, { £f,(p-" + p,p^) } 

+ ^J2^'-{ ^M«(P™- + P^P^ + pITp,) } + (p ^ ^) • (10) 
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As one would expect, this equation contains the three-point correlator p™", and similarly 
the evolution equation for p'^^^J^x contains the four-point correlator etc. Consequently, 
one can never exactly solve this set of equations, truncated at any finite order. 

However, the limit Z ^ 1 facilitates an approximate solution: Let us imagine that 
we start from an initial state p™ = ^^P^n without any correlations (i.e., p'j^J^ = and 
P/J°K = O5 6tc.) such as the Mott state ([3]). In this case, the right-hand side of fITOl) 
is suppressed by 0{1/Z) and thus the time evolution creates only small correlations 
idtp'T". Moreover, if these correlations are small initially p'^^" = 0{1/Z), they remain 
small - at least for a finite amount of time - because there is no term in (TTOj) to 
compensate the 0{1/Z) suppression. Note that the sum ^^ in ( ITOl) might scale with 
Z, but this is compensated by the 1/Z factor in front of it. On the other hand, if we 
insert p"^"""" = 0{1/Z) into Q), we find that we can neglect this term and arrive at an 
approximate equation containing on-site density matrices only 

^^^^'^K{ClJf,p^]+ C^p^, (11) 

The approximate solution fp of this self-consistent equation is valid to lowest order in 
1/Z, i.e., p^ = p°-|- 0{1/Z) and reproduces the well-known Gutzwiller ansatz [1811211135] . 
If we now insert this approximate solution fP into f lTO|) , we get an approximate evolution 
equation for the two-point correlator 

^dtp^- = c,p-' + 1 c,^pIpI + ^J2^'^{ ^UpITpI + pITpI)} 

- f tr, { C^jy^} + {p^u)+ 0{1/Z') . (12) 

Note that we have assumed that the three-point correlations p™''^ do not spoil this line of 
arguments and are suppressed by 0{1/Z'^) in complete analogy. This is indeed correct 
and can be shown in basically the same way, see Appendix [171 More generally, we find 
that £-point correlations are suppressed as 0{1/Z^~^), i.e., 

P, =0 (Z°) 

P-" =0{1/Z) 

(f^ = o {l/Z') 

p;°"a = O [l/Z^) , (13) 

and so on, see Appendix [T71 This hierarchy (fT3|) is related to the quantum de Finetti 
theorem [36], the generalised cumulant expansion [37], and the Bogoliubov-Born-Green- 
Kirkwood-Yvon (BBGKY) hierarchy [32], but we are considering lattice sites instead of 
particles. As an example for the four-point correlator, let us consider observables A^, 
Biy, Cki and Dx at four different lattice sites, which have vanishing on-site expectation 
values (A^) = {13^) = {C^) = (Dx) = 0. In this case, the hierarchy ( 1T3|) implies 

{A^B^C^Dx) = {A,A){C.Dx) + {A^C^){B,Dx) + {A,Dx){B,C^) 



Correlations in the Bose & Fermi Hubbard Model 7 

+ O [l/Z') , (14) 

which resembles the Wick theorem in free quantum field theory (even though the 
quantum system considered here is strongly interacting). 

4. Mott Insulator State 

Now let us apply the hierarchy discussed above to the Bose- Hubbard model ([1]). To 
this end, we start with the factorising Mott state (j3]) at zero hopping rate J = as our 
initial state 



r = Q^p'^ = Q^\i),{i\- (15) 

Then we slowly switch on the hopping rate J{t) until we reach its final value. In view 
of the finite energy gap, the adiabatic theorem implies that we stay very close to the 
real ground state of the system if we do this slowly enough. Of course, we cannot cross 
the phase transition in this way (i.e., adiabatically) since the energy gap vanishes at the 
critical point, see Section |5] below. 

Since we have (6^) = in the Mott state, Eq. ( ITTl) simplifies to 

^dtP, ^^Y^^^A ^lJ,P^] + ^,h = -- pj = |1),,(1| . (16) 

Thus, to zeroth order in 1/Z (i.e., on the Gutzwiller mean-field level), the Mott insulator 
state fp for finite J has the same form as for J = 0. To obtain the first order in 1/Z, 
we insert this result into ( IT2l) . Again using {b^) = 0, we find 

+ ^ E ^^'^ { ^'«prp° + ^Lpt^'pI] + o{i/z') . (17) 

Formally, this is an evolution equation for an infinite dimensional matrix ff°". 
Fortunately, however, it suffices to consider a few elements only. If we introduce 
p^ = |1)„(2| and h^ = |0) (1| as local particle and hole operator^, all the interesting 
physics can be captured by their correlation functions (for fi ^ u) 



C = HK) = tr {phlK] = tr,. {p'.TKh 
III = (Kp^) = tr {phlpl = tr,. {p^r/^iA 



ful = (PJ./^.) = tr IppIk] = tV {pTPih 



fil = iPlPu) = tr {pPlPu} = tr,. [pITpIPu] ■ 
To first order in 1/Z, these correlation functions form a closed set of equations 

J 



I These excitations are sometimes [inilMlllO] called doublons and holons. 
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.12 JV2^ 



+ Uflt - -^T,, (19) 

J 



^dtC = + ^ E (^-(/S + ^O + ^T,M^' + ^i 



12^ 
KU/ 



Ufil + "-^T,, (20) 



.21 , ■J'^r 



^dJll = ^^,C = -^ Yl iT,.C - T.Jil) . (21) 

This truncation is due to the fact that the correlation functions /™" involving higher 
occupation numbers m > 3 or n > 3 do not have any source terms of order 1/Z 
and hence do not contribute at that level. Exploiting translational symmetry, we may 
simplify these equations by a spatial Fourier transformation with 

T,u = ^Y.'^^'^'^'^''''^^^ (22) 

k 

tf = ^E/k'e^'■^''''^'"'^ (23) 

k 

where A^ denotes the number of lattice sites (which equals the number of particles in 
our case). Formally, in order to Fourier transform equations ( iT9ll2Tl) . one should add 
the summands corresponding to n = fi and k, = u. Since these terms are of order l/Z"^, 
they do not spoil our first-order analysis. However, when going to second order 1/Z'^ 
(see Section [9] below), they have to be taken into account. 

After the Fourier transformation (!22|) and (123|) . Eqs. (I19H21I) become 

(,5^ _ f/ + 3JT^)Jl' = - V2JUfi' + /f + 1) , (24) 

{^^t + U- 3JTk)/f = + V2JTk(/^i + /f + 1) , (25) 

tdJl' = id,}f = V2JT^UV - ff) ■ (26) 

From the last equation, we may infer an effective particle-hole symmetry /^^ = /^^ 

valid to first order in 1/Z. With this symmetry, any stationary state (such as the 

ground state) with dtf^^ = must obey the condition 

12 _ ,21 _ V2JU2f^' + 1) 
/k -/k - U-3JT^ ■ ^ '' 

The remaining unknown quantity f^ can be obtained in the following way: The 
evolution equations ( 124H26|) leave the following bilinear quantity invariant 

5t[/k'(/k' + l)-/kX]=0' (28) 

which remains valid even for time-dependent J{t). Thus, starting in the Mott state 
on]) at zero hopping rate J = with vanishing correlations f^{t = 0) = 0, we get the 
additional condition f^{f^ + 1) = f-^f^ for all times t > 0. Thus, to first order in 
1/Z, the ground state correlations read (for ^ ^ v) 

\'^^"'v/ ground = \P/^Pi/ / ground = "T7 / ^ T~ 6 '' " (29) 

k 
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(/^iP.)ground = {piK),ronna = ^ ^ ^^C^^'^— ) . (30) 



Here we have used the abbreviation 



k 



c^k = yf/^ - 6Jf/Tk + .PT^ , (31) 

which is just the (non-zero) eigenfrequency of the evolution equations (I24ti26p for non- 
stationary states and will become important in the next Section. 

The above equations (129|) and (1301) describe the correlations and are valid for fi ^ i' 
only. The correct on-site density matrix p^ can be obtained from ([9]) which shows that 
non- vanishing correlations lead to small deviations from the lowest-order result p^ . As 
one would expect, the quantum ground-state fluctuations manifest themselves in a small 
depletion of the unit-filling state p*^ = |1)„(1| given by a small but finite probability for 
a particle /2 = tr{p^|2)^(2|} or a hole /o = tr{p^|0)^(0|}. To first order in 1/Z, we get 
from (ED 

^^Jo = ^^J. = Y: ^^UV - fl') = 4 E ^*/^^ ' (32) 

k k 

where we used equation (l26l) in the last step. This equation can be integrated easily 
and with the initial conditions /o(t = 0) = /2(t = 0) = we find the 1/Z-corrections to 
the one-site density matrix 

(p%) = (KK) =h=k=^j: "~ '2 J" " "" ■ (33) 

k ^ 

Note that, even though the right-hand side of the above equation looks like that of (12^ 
for fi = u, one should be careful as they are derived from two different equations: iQ 
and (im]). 



5. Mott— Superfiuid Quench 

After having studied the ground state properties of the Mott phase, let us consider a 
quantum quench from the Mott state to the superfiuid regime. This requires a time- 
dependent solution of the evolution equations ( I24ll26l) which crucially depends on the 
eigenfrequency (!3T]) . In view of the definition (122|) . T^ adopts its maximum value 
2k=o = 1 at k = 0. Thus aJk=o = ^^ corresponds to the energy gap of the Mott 
state mentioned in Section O For J = 0, we have a fiat dispersion relation u^ = U. 
If we increase J, the dispersion relation Wk bends down and the minimum at k = 
approaches the axis. Finally, at a critical value of the hopping rate [12] 

Jcrit = f/(3 - VS) , (34) 

the minimum ci;k=o touches the axis and thus the energy gap vanishes AS = 0. This 
marks the transition to the superfiuid regime and we cannot analytically or adiabatically 
continue beyond this point. However, nothing stops us from suddenly switching J to 
a final value Jout > <^crit beyond this point. Of course, this would not be adiabatic 
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10 



anymore and we would no longer be close to the ground state. For hopping rates J 
which are a bit larger than the critical value J > Jcnt, the dispersion relation dives 
below the axis and the u^ become negative for small k. Thus, the eigenfrequencies Wk 
become imaginary indicating an exponential growth of these modes, i.e., an instability. 
This is very natural since the quantum system "feels" that the Mott state is no longer 
the correct ground state. 

If we consider even larger J, we find that the original minimum of the dispersion 
relation w^ at k = splits into degenerate minima at finite values of k when J = 3U, 
while k = becomes a local maximum. This local maximum even emerges oj^^q > 
on the positive side again for J > f/(3 + v^)- Nevertheless, there are always unstable 
modes for some values of k, see Fig. [1] and compare 
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Figure 1. Dispersion relation w^/t/^ in one dimension for different values of J/U. 



After these preliminaries, let us study a quantum quench from the Mott state to 
the superfluid phase. For simplicity, we consider a sudden change of J(t) = JQ(t) from 
J = to the final value of J (but the calculation can easily be generalised to other 
scenarios). Solving the evolution equations ( l2Hl26l) for this case, we find 

— V 4 J^T^ - ~ ^°^^^kt) „ik.i 



\ '''/J '''!^/ quench \P;^Pi// quench 



N 



■'\^l-i ^u ) 



uJt 



(35) 



{hlPXnenc, = ^ ^ V2JUU - 3JT^) ^ ^"^M) e^k- 



(Xp-Xy) 



■ (.X^ Xjyj 



(36) 



The remaining correlation can simply be obtained via {plh ) = ((/itp^))*- The correlator 
in terms of the original creation and annihilation operators 6j, and b^ is just a linear 
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combination of these correlation functions 

k ^ 

Note that the momentum distribution 

^W = ]^E^*'''^''*""^(^i^-) (38) 

fJ,U 

which is basically the Fourier transform of {b'^bv)-, can be measured by time-of- flight 
experiments [HI SI] The quench J{t) can be realised experimentally by decreasing the 
intensity of the laser field generating the optical lattice (which lowers the potential 
barrier for tunnelling and thus increases J). Thus the above prediction should be testable 
in experiments. 

As explained above, after a quench to the superfiuid regime, the frequencies cuk 
become imaginary for some k and thus these modes grow exponentially. As a result, the 
expectation value will quickly be dominated by these fast growing modes and so most 
of the details of the initial state will become unimportant. Of course, this exponential 
growth cannot continue forever - after some time, the 1/Z-expansion breaks down since 
the quantum fluctuation are too strong and the growth will saturate. 

6. Equilibration versus Thermalisation 

Instead of a quench from the Mott to the superfiuid phase, we can also study a quench 
within the Mott regime. Again, we consider a sudden change of J from zero its final 
value for simplicity - but now the final value J lies below the critical point J < Jcrit- In 
this case, all frequencies are real cuk G M and thus there is no exponential growth - all 
modes oscillate. Apart from this point, we can use the same solution as in ( I35ll37|) . For 
an infinite (or at least extremely large) lattice, the oscillations in ( I35ll37|) average out for 
sufficiently large times t and thus these observables approach a quasi-equilibrium value 

1 A J2rn2 

{Kh.U^^ = (p!.P.)cqun = ^J2 —2^ e^^-(— ) (39) 

^ k ^^ 

The quasi-equilibrium values for the local (on-site) particle or hole probability can be 
derived in complete analogy to the previous case 

1 4J^T^ 

(Pip^)cqmi = {h^hl)^q^ii = 1^Y1 ~7l2^ ■ (^^) 



k 



k 



Having found that the observables considered above approach a quasi-equilibrium state, 
it is natural to ask the question of thermalisation. This is one of the major unsolved 
questions (or rather a set of questions) in quantum many-body theory ^5[ B6l W7\ HH HO] 
In one version, this question can be posed in the following way: Given an interacting 
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Figure 2. Time-dependence of the depletion {p\p^) and the nearest-neighbor 
correlations functions {n^Jiu) and {p\Pu) in three dimensions after the quench within 
the Mott phase from J/U = to J/U = 0.14 in comparison to their ground state 
values. After quasi-equilibration, (p]^p^)quench and (pj^Pi.) ground as well as (pJjP^)quench 
and (ptp/i)ground differ roughly by a factor of two. 



quantum many-body system (for example the Bose-Hubbard model) on an infinite lattice 
in a appropriately excited state (such as after a quench), do all localised observables 
(e.g., {pjj_p^) and (h'^^h^,)) settle down to a value which is consistent with a thermal state 
described by a suitable temperature? 

Even though we cannot settle this question here, we can compare the quasi- 
equilibrium values obtained above with a thermal state. To this end, we derive the 
thermal density matrix p^ corresponding to a given (inverse) temperature /cb^ = 1//^- 
Using the grand canonical ensemble, the thermal density operator is given by 

P^ = tr{e-^iH-,m} ' ^^'^ 

where chemical potential /i will be chosen such that the filling is equal to unity. Since 
we cannot calculate p/j exactly, we employ strong-coupling perturbation theory, i.e., an 
expansion in powers of J. It is useful to introduce the operator [51] 
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where Hq is the diagonal on-site part of the grand canonical Hamiltonian H — fiN and 
Hi is the hopping term. This operator satisfies the diff^erential equation 

d0R{/3) = -Hi{t3)R{(3) , (44) 

where Hi{f3) = e^^^Hie'^^'^. In analogy to time-dependent perturbation theory, the 
operator R can be calculated perturbatively by integrating this equation with respect 
to /3. In first-order perturbation expansion, we have 

''' - ^^IF^} (^ ^ 1 5 r.« u(,^ - .. - 1) "'•"" ^ ^(^')) ■ («) 

Obviously, the correction to first order in J does not affect the on-site density matrix 
p^ but the two-point correlations. Thus, we find that the quasi-equilibrium state of the 
on-site density matrix p^ can indeed be described by a thermal state provided that we 
choose the chemical potential as /i = U/2 which gives 

P-I3U/2 p-fiU/2 

pM ^ —^ I0),(0| + (1 - e-^^/2) |1)^(1| + ^- |2)^(2| . (46) 

The particular value /x = U/2 of the chemical potential ensures that (in first order 
thermal perturbation theory) we have on average one particle per lattice site and the 
particle-hole symmetry {f)\iP^i) = (^t^^t)- To obtain the correct probabilities, we have 
to select the temperature according to 

k ^ 

which can be deduced from ( 14T|) and (146|) . Since the depletion is small {Pl^p^) = 0{1/Z), 
we obtain a low effective temperature which scales as T = 0{U/ In Z). Accordingly, 
consistent with our 1/Z-expansion, we can neglect higher Boltzmann factors such as 

Of course, the fact that the on-site density matrix p^ can be described (within our 
limits of accuracy) by a thermal state does not imply that the same is true for the 
correlations. To study this point, let us calculate the thermal two-point correlator from 
(gS]). To first order in J and 1/Z = 0{e-^^/2), we find 

ChlP^)^ = {pIK)p = ^^^^ + 0{J') + 0{llZ') , (48) 

while {li^^h^)p and {pltPu)i3 vanish (to first order in J). If we compare this to the quasi- 
equilibrium value {hjj^Pi,)cqui\ in fHOj) . we find that they coincide to first order in J 

{hiPuUun = {piK)e,u, = ^^^^ + 0{J') + 0{l/Z') . (49) 

This is perhaps not too surprising since the same value can be obtained from the ground 
state fiuctuations {h'^^Py) ground = {P^hy) ground in (j30|) after expanding them to first order 
in J. Due to the low effective temperature T = (9(f//lnZ), the lowest Boltzmann 
factor is suppressed by e'^^/"^ = 0{1/Z). As a consequence, because the correlations 



2(pip,)equn = ^Y1 ^^ = ^(V^) , (47) 
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are small 0{1/Z), their finite-temperature corrections are even smaller 0{1/Z^), and 
thus can be neglected. 

The same is true for the other correlations {h^Ji^) = (pip^)- All of them: the ground 
state correlators (^|j/ij,)ground = {p^^Pu) ground in ( 129|) . the quasi-equihbrium correlators 
{%h^,)cqmi = {pli.Pu)eqmi in ( l39|) . as wcU as the thermal correlators {h)Jii,)ii and {p^^Pi) p 
vanish to first order in J. Therefore, to first order in J and 1/Z, the thermal state can 
describe the observable under consideration. However, going to the next order in J, this 
description breaks down. This failure can even be shown without explicitly calculating 
i?(/9) up to second order. If we compare the quasi-equilibrium correlators in ( l39l) 



4 P 

{hiK)equil = (p;A)equil = JJ^ Y. ^^^^^^ + ^^^'^ + ^^ V^') , (50) 

with the ground state correlations in (!29|) . expanded to the same order in J 

2/2 

(/^i/^.)ground = (]5i]5jg,o,„d = JJ^^Y.^,nT^. + 0{j') + 0{l/Z') , (51) 

K 

we find a discrepancy by a factor of two. I.e., after the quench, these correlations 
settle down to a value which is twice as large as in the ground state. This factor of 
two has already been found elsewhere in the context of standard time-dependent and 
time-independent perturbation theory, see also [50] . This is incompatible with the small 
Boltzmann factors e'^^/"^ = 0{l/Z) and would require a comparably large effective 
temperature T = 0{U) instead of T = 0{U/ In Z). However, such a large effective 
temperature T = 0{U) is inconsistent with the small on-site depletion P7|) . 

In summary, the considered observables settle down to a quasi-equilibrium state 

- but this state is not thermal. Thus, real thermalisation - if it occurs at all - 
requires much longer times scales. This seems to be a generic feature and has been 
discussed for bosonic [HH |52l [53] and fermionic systems [HH [551 EEl [571 EHl [59] and 
is sometimes called "pre-thermalisation" . This phenomenon can be visualised via the 
following intuitive picture: The excited state generated by the quench can be viewed 
as a highly coherent superposition of correlated quasi-particles. During the subsequent 
quantum evolution, these quasi-particles disperse and randomise their relative phases 

- which results in a quasi-stationary state. However, the quasi-particles still retain 
their initial spectrum (in energy and quasi-momentum) , which could be approximately 
described by a generalised Gibbs ensemble (i.e., a momentum- dependent temperature). 
In this picture, thermalisation requires the exchange of energy and momentum between 
these quasi-particles due to multiple collisions, which changes the one-particle spectrum 
and takes much longer. Ergo, one would expect a separation of time scales - i.e., 
first (quasi) equilibration and only much later thermalisation - for many systems in 
condensed matter, where the above quasi-particle picture applies. 
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In the following, we study the impact of a spatially constant but possibly time- dependent 
force on the particles, which could correspond to a tilt of the lattice, for example 
[60l|6ll[62l|63l[Mll65l|66l|67l|68l[Ml[7Ql[7ll[^ This scenario can be 

described by a generalisation of the Hamiltonian ([T]) 



H 



J 



U 






n^[n^. 



+Ef. 



fi^fj^ . 



(52) 



The external potential V^(t) = x^ ■ E(t) will be identified as an effective electric field 
E(t) and will be time-dependent in general. If we insert this modified Hamiltonian into 
f lT6|) . we find that the potential V^ has no effect to zeroth order 0{Z^), i.e., the solution 
Pu = |1)„(1| remains the same. In other words, the Gutzwiller mean field is not affected 
by the tilt in the Mott state (in the superfluid phase, this would be different). However, 
the next-order 0{1/Z) quantum correlations can lead to the creation of particle-hole 
pairs via tunnelling over one or more lattice sites. 

In order to study this effect, let us generalise the evolution equations flTOl) and fl2Tl) 
in the presence of the potential V^ 



{tdt + v^-v,-u) f 



12 



/IK 



{id, + v,- K) fll 



z ^ ' 



'^fll + ^Vll + V2C 



flK 



/ .21 _ H2N 



(53) 
(54) 



and the same for f'^l, such that we again have an effective particle- hole operator 
symmetry fj^l = /^^ to lowest order in 1/Z. Here we have already used translational 
invariance. The tunnelling probability can now be obtained by solving the above 
equations. However, instead of solving them directly, we can simplify the problem by 
effectively factorising these equations: If we introduce the effective differential equations 
for hf, and p^. 



idt - 


--t] 


idt- 


u' 



Pf^ 



K 



-Tt 



■Pu 



-Tt 



- hv 
2 



+ V2K 

V2pu 



(55) 



and exploit the initial conditions {W^hy)^ = 5^j, and {h)^'Pv)o = (pJj^i/)o = {p\Pu)o = 
valid in the Mott state, we exactly recover (153|) and (l54l) to first order in 1/Z. 

For potentials of the form V^(t) = x^ ■ E(t) it is possible to apply the Peierls 
transformation and absorb the potential in a phase. After the Fourier transformations 



h^{t) = exp <^ -i / dt'Vf,{t') \ ^ h^{t) exp{ik ■ x^} 



(56) 
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p^it) = expl-t I dt%{t')\j2pUt) exp{^k ■ x^ , (57) 

r^.(t) = ^5Z^k(t)exp|2k-(x^-x,)+^ /" dt'[V^it')-V,it')]\ ,(58) 

the effective evolution equations ( l55l) become 

tdtk = +1 [3 JTk(t) -U]h^ + v^JTk(t)]3k , (59) 

idtPk = -I [3Jrk(t) - f/] Pk - V2 JTk(t)/zk . (60) 

Note that Tk(t) exphcitly depends on time here and this time- dependence encodes the 
potential V^(t). Introducing the effective vector potential A(t) which generates the 
effective electric field E(t) in V"^(t) = x^ ■ E(t) via E(t) = dtA(t), this is equivalent to 
the substitution k — )■ k + A(t) in complete analogy to the minimal coupling procedure 
Tk(t) = Tk+A(i) known from electrodynamics. 

The most general solution of ( l59l) and ( l60l) can be written as 

/ik = /k+(t)ik + /kW5k, (61) 

p^=g^{t)A^ + g^{t)B^, (62) 

where Ay^ and i?k are time-independent operators while f^ and g^ are time-dependent 
c-number functions. In analogy to the previous case, we assume that we start in the 
Mott state ([3]) with J = V^ = 0. In this case, the equations (159|1 and (l60l) decouple and 
we may choose Ak = ^k ^^'^ -^k = Pk which imply {Al_Ap)o = 5k,p and {BI.Bp)q = 0, 
as well as, f^{t) = exp{iUt/2) and g^{t) = exp{—iUt/2) while the other two vanish. 

Now we imagine the following sequence: First we switch on J adiabatically, then 
we apply the potential V"^(t) for a finite period of time, and finally we switch off J 
adiabatically. Thus, at the very end, the equations ( 159|) and (!60l) decouple again and 
the final particle operator p^^ oscillates with positive frequencies exp(— if/t/2) while 
the final hole operator h^^ oscillates with negative frequencies exp{+iUt/2). However, 
during the time-evolution, positive and negative frequencies will get mixed in general 
by the time-dependence of Tk(t) = Tk+A(i), i-e., the potential V^(t). Thus, the initial 
and final particle/hole-operators will be connected by a Bogoliubov transformation 

^r = «kPL"+/3k/^L", 

hr = aX + /3kPL" , (63) 

where the Bogoliubov coefficients ak and /3k satisfy the relation |akP — |/3kP = 1- In 
view of the initial conditions (A[.Ap)o = 5k,p and (5^5p)o = 0, we find 

(pLPk)out = |/3kr, (64) 

which gives the probability to create a particle in the mode k. Since particles (i.e., 
doubly occupied lattice sites) and holes (i.e., empty lattice sites) are always created in 
pairs, we get the same probability for the holes. Note that k is the canonical and not 
necessarily the mechanical momentum due to the substitution k — )■ k + A(t) mentioned 
above. 
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8. Analogue of Sauter-Schwinger Tunnelling 

The precise amount of mixing which determines the Bogoliubov coefficients a^ and /3k 
can be derived from the evolution equations fl59|) and fl60|) . Turning these two first-order 
differential equations into one second-order equation, we get for ^fj^ and f^, 

d^fi-^dJ^+\'^ + ^^]f^ = 0, (65) 

dh^ - ^i^d^gi +{^-^^\9^ =0. (66) 




4 2Ti 



k 

With the substitutions f^ = -\/7kUk and g^ = y/T^Vk, we may eliminate the ffist-order 
terms and arrive at 

^2 fiol Uf^ fk 3t2 , 



cot .UT^ Tk 'JTl 
2Tk 2Tk 4r,? 



52„, , ^k .^-'k , Jk ^Jk 



^X+|^-^W^ + ;^-^h'k=0. (67) 



Now we consider a small tilt of the lattice, corresponding to a weak electric field |E| ^ f/. 
In this case, we may approximate the above equations by neglecting the terms with Tk 
and Tk since Tk = 0(E). Furthermore, for a weak tilt, the particles have to tunnel 
across many lattice sites in order to gain enough energy and to be able to overcome 
the energy gap before a particle-hole pair can be created. Thus, we may consider large 
length scales, corresponding to small wavenumbers k and Taylor expand the Tk(t) 

Tk(t) = Tk+AW = 1 - e[k + A(t)]2 + (9(k4) , (68) 

where ^ is the stiffness. With these approximations, we find that (!67l) simplify to 

^U^. + (me\4 + 4[k + A{t)f) 0k = , (69) 

which is just the Klein- Fock-Gordon equation describing charged scalar particles in an 
external electromagnetic field, provided that we identify the effective speed of light 

4 = \ A^u - J) , (70) 

while the effective mass is given by half the energy gap AS 

ml^ct^ = ^iU'-QJU + j'). (71) 

Consequently, there is a quantitative analogy between the tilted Bose-Hubbard lattice 
and the Sauter-Schwinger effect, i.e., electron-positron pair creation out of the quantum 
vacuum due to an external electric field, sketched in Fig. [3] and the following table: 
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Sauter-Schwinger effect 


Bose-Hubbard model 


electrons & positrons 


particles & holes 


Dirac sea 


Mott state 


mass of electron/positron 


energy gap AS 


electric field E 


lattice tilt V^ 


speed of light c 


velocity c^s 



a) 



+mc -o - ^- - o 
^2 . . . 



-mc 



£ 

-O O O - 



c) 




b) 



j»-<' 




Figure 3. Sketch of the analogy: a) Dirac sea for E ^ 0, h) Sauter-Schwinger 
tunneling for E y^ 0, c) Mott state with energy gap A£, d) tunneling in tilted lattice. 

We can now use this analogy to apply our knowledge of the Sauter-Schwinger effect 
[76l [771 [78l [79l l80| ISTj to particle-hole creation in the tilted Bose-Hubbard model ~ 
as Richard Feynman said: The same equations have the same solutions. For example, 
consider a purely time-dependent electric field of the following form 



E 



cosh^ (t/r) 



(72) 



Such a profile is called Sauter pulse since F. Sauter was the first to realise (already in 
1931) that the Dirac equation and the Klein-Fock-Gordon equation in the presence of 
such a field can be solved exactly in terms of hypergeometric functions (although he 
considered the form with t and x interchanged). From the exact solution of the scalar 
field case, one obtains [821 [83] 

cosh (7rr[a;+ — u^]) + cosh (^r^/AE^c 



n2^2^4 



l/3k| 

with the abbreviations 

u± - 



2 sinh(7rra;+) sinh{7rTUJ- 



(73) 



c^ihTEory + mlc'^ + klc^. (74) 

Here k± denotes the momentum perpendicular to the electric field and m^ is the mass of 
the electron. Via the analogy established above, expression ( !73|) yields the momentum 
dependent particle-hole pair creation probability in a Mott state subject to a time- 
dependent tilt according to Eq. (172|) . For various pulse lengths r, this result plotted in 
Fig. [Hland compared to numerical results for a one-dimensional Bose-Hubbard lattice. 
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In the static limit r — ;■ oo, Eq. ([72]) reproduces the well-known expression 



|/3k 



exp 



-7r- 



Eoc 



(75) 



As we see, the electron-positron pair creation probability is exponentially suppressed for 
small electric fields Eq. Inserting the translation formula (ITOll and (ITTll . we get the same 
exponential suppression for the particle-hole pair creation probability via tunnelling in 
tilted Mott lattices. Thus, in order to actually verify this prediction experimentally, the 
tilt should not be too small. In this case, the terms with Tk and Tk we have neglected 
earlier due to Tk = 0(E) might start to play a role. Thus, let us estimate the impact 
of these contributions. Including the terms involving Tk and Tk, we find 

dtU], + [mlgcls + klclg + clsik, - Eotf] u^ 

+^ [-E^ + iEoU{k, - Eot)] Uk = , (76) 

where we have assumed a constant field (r — ?■ oo) for simplicity. The above differential 
equation can be solved in terms of parabolic cylindrical functions from which the pair 
creation probability is determined to be 



|/3k 



exp 



n 



EqC, 



cff 



2 4 



+ clskl - ^El + e 



2 EIU' 



(77) 



In Fig. m we depicted the dependence of the particle- hole creation probabihty (pjj5^) 
Xlk l/^kP/A^ on the potential gradient. 




Figure 4. Dependence of —\\\{{p^p^)) on the lattice tilt for J/U = 0.1. The black 
line represents the standard result ([75]) for the static Sauter-Schwinger effect while the 
red curve deviates due to perturbative corrections in Eq given by the lattice structure, 
see equation fTT]). 



9. Second Order in 1/Z 

So far, we have only considered the first order in 1/Z. Now let us discuss the effect of 
higher orders by means of few examples. Let us go back to the derivation from ( ITO|) 
to (fT2|) and include 1/Z^ corrections. To achieve this level of accuracy, we should not 
replace the exact on-site density matrix p^ by it lowest-order approximation p^ but 
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include its first-order corrections in ( 15^ . i.e., the quantum depletion /o = 0{1/Z) of 
the unit filling (Mott) state. This results in a renormalisation of the eigenfrequency 



<" = \JU' - 6JTk(l - 3/0) + J^Tlil - 3/0)2 , (78) 

which indicates a shift of the Mott-superfiuid transition to slightly higher values of J, 

J^Il = U ^f^ > Ui3 - 2^2) , (79) 

see Appendix [191 Since the net effect can roughly be understood as a reduction of the 
effective hopping rate J^™ = J(l — 3/o), it is easy to visualise that this implies also a 
decrease of the effective propagation velocity. 

There are also other l/Z"^ corrections in flT2l) such as the three-point correlator p"^" 
but they act as source terms and do not affect the eigenfrequency (at second order). 
However, there are other quantities where these source terms are crucial. In particular 
we consider correlation functions which are of the form 

Faifi,iy) = {d,d,)-{d,){d,), (80) 

and vanish to first order in 1/Z, in contrast to the off-diagonal long-range order 
(aj^Oiy) discussed above. One important example is the particle-number correlation, 
i.e., {n^fijj) — 1. After a somewhat lengthy calculation, we find for the ground state 
correlations (see Appendix [T9l) 

^^(/^' -) = ]| E e^^-^->^-^-'^^ [fl'f'^ - fl'fl') , (81) 

p.q 
where f^, f"^ and f^ are given through equations fl271) and fl28l) . It is also possible to 
calculate this quantity via a perturbation expansion into powers of J /U . Note, however, 
that the above result is not perturbative in J /U , see, for example, the non-polynomial 
dependence of cjk on J. 

These predictions could be tested experimentally by site- resolved imaging, i.e., 
measurements on single lattice sites [101 IHH ESI ES]. In some of these experiments, 
the particle number per site is not directly measured, but only the parity - i.e., whether 
the number of particles on a given lattice site is even or odd [1^ . The parity correlator 
reads 

i^(-i)"(/^, -) = ^.Y. e^^-^->^-^-^^ {fl'fll + /;Vf ) . (82) 

p,q 

Again, this result can be compared with a perturbative expansion into powers of J /U . 

Assuming a hyper-cubic lattice Z^ in D dimensions with nearest-neighbour hopping 

(i.e., Z = 2D), we obtain up to quartic order in J 

F(_i)n(/i,z/)= (-^j 8n{n + l) + 

+4 - 22Z + 9^2] + C( J6) ^ (83) 
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where /i and z/ are nearest neighbours and n = (n^) is an arbitrary (integer) filhng. 
Inserting n = 1 and keeping only the lowest-order l/Z"^ terms, we may compare this 
result with ( l82l) . after an expansion into powers of J, and find perfect agreement. 

However, there is an interesting observation regarding the above equation: In one 
spatial dimension with Z = 2 nearest neighbours, the J^ contribution in the above 
equation is negative. This suggests that the parity correlator assumes its maximum at 
a finite value of J (in the Mott phase), which can indeed be confirmed by numerical 
simulations, see, e.g., [20] and Section [101 In two or more spatial dimensions, the 
situation is different. Even though the parity correlator should still assume its maximum 
at some finite value of J, this value is quite close to the phase transition or already in 
the superfluid regime. Thus, this maximum is not visible in our 1/Z-expansion starting 
in the Mott state, which predicts a monotonously increasing parity correlation in its 
region of applicability. 

In analogy to Sections [5] and [6l we can also study the correlations after a quantum 
quench with J{t) = JQ{t). Again, there are no contributions to the particle-number 
and parity correlations in first order 1/Z - but, to second order 1/Z, we find formally 
the same expressions as in the static case ( l82l) and ( IHTI) 



^"(/^'^) = ^$^e'(P+^Hx.-x.) (/n(i)jn(^) _ /^^(t)/^ (t)) , (84) 



Ar2 
p>q 

and 

Fi-,Af^,u) = A^e^(P+q)-(x.-.) {fl\t)f^\t) + fl\t)fl\t)) , (85) 
p,q 
where f^{t)J'^{t) and f^{t) are now given by equations fl5^ and fl5Bl) . The parity 
correlations after a quench have been experimentally observed in a one- dimensional 
setup [in]. Although the hierarchical expansion relies on a large coordination number, 
we find qualitative agreement between the theoretical prediction fl85|) for Z = 2 and the 
results from [TH]. For large times t and distances x^ — x^, we may estimate the integrals 
over p and q in the expressions (!84l) and (185|) via the stationary-phase or saddle-point 
approximation. The dominant contributions stem from the momenta satisfying the 
saddle-point condition 

Vk [k ■ (x^ - x,) ± u^t] = . (86) 

Thus their structure is determined by the group velocity Vk = Vk^k- If the equation 
x^ — x^ = ±Vkt has a real solution k, i.e., if the distance x^ — x^ can be covered 
in the time t with the group velocity Vk, then we get a stationary-phase solution - 
otherwise the integral will be exponentially suppressed (i.e., the saddle point k becomes 
complex). For small values of J, the maximum group velocity is given by vj^*^^ ^ 3J, 
which determines the maximum propagation speed of the correlations, i.e., the effective 
light-cone. 
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10. Numerical Simulations for the Bose-Hubbard model 

In the following we analyze the Bose-Hubbard system ([T]) numerically in one and two 
dimensions. All calculations are carried out on a finite lattice with L lattice sites and 
N bosons. 

10.1. General formalism for the one- dimensional Hubbard model 

The eigenstates of lattice systems are calculated by means of exact numerical 
diagonalisation of the Hamiltonian matrix [871[HHl[H9ll9nil9ll[92l|93l|9l[95] which 
can be obtained from the Hamiltonian ([1]) using the basis of Fock states 

|nr) = (g) |nr,) , r = l,...,V, V = ^^,|/__~ |f , (87) 

where F labels the configuration of the bosons and the occupation numbers of individual 
lattice sites n-p^ satisfy the condition 

L 

N = J2^rf^- (88) 

The matrix dimension can be reduced by factor L for homogeneous lattices with periodic 
boundary conditions (fei+i = 6i). In this case the Hamiltonian commutes with the 
unitary translation operator T which acts through cyclic permutation on the lattice 
bosons. Due to the periodic boundary conditions, the operator satisfies T^ = 1. As a 
basis one can use linear combinations of the Fock states flHTD in the form 



\^Kr)=AfrJ2[-^) l^r), (89) 

which are eigenstates of the operator T for the eigenvalue tk = exp (iK) . A/r are 
normalisation constants chosen such that (nxr|nx'r') = ^rr' ^kk'- Here the state |nr) 
cannot be obtained by cyclic permutations of |nr') with F' ^ F. The eigenstates of the 
Hamiltonian have the following form 

Vk L-1 

\Kn) = J2CKnr\nKr), n = l,...,VK, Y.^^ = ^^ (90) 

r=i K=o 

and the corresponding eigenenergies are denoted by Exn- 

If the complete eigenvalue problem is solved, one can work out expectation value 

of any operator O at arbitrary temperature in a canonical ensemble as 

(^)^ = Ym E^^^^^l^l^^^^) ^^p (-^) (91) 

with the partition function 

Z(r) = $:exp(-:^). (92) 
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Using the set of the basis states (I89p we were able to solve numerically the complete 
eigenvalue problem for N = L = 9. 

In order to study zero-temperature properties of the system, it is sufficient to 
calculate the ground state. This can be done exactly with the aid of iterative numerical 
solvers for sparse matrices of large dimensions along the lines of [91] . We were able to 
do this up to iV = L = 14. 

10.2. Energy spectrum 

In the limit of vanishing hopping, J = 0, the basis states ( l89l) are the eigenstates of the 
Hamiltonian ([1]) which are, apart from the ground state, degenerate. The ground state 
has equal occupation numbers at each lattice site, that is nr^ = n. It exists only at 
K = and has the energy 

Eoi = Ljn{n - 1) . (93) 

The energy eigenvalues of the degenerated excited states are given through 

L TJ 

^^r = $^lTnrMKM-l) (94) 

and do not depend on K, corresponding to flat energy bands. The lowest band contains 
L{L — 1) degenerate eigenstates with the energies Ekv = -E'oi + U. These states 
correspond to bosonic configurations with the same occupation numbers n at any site 
except two, one of which contains n — 1 boson and the other one n + 1. The highest 
band contains L degenerate states with all atoms sitting at one lattice site. These states 
have the energy Ekv = UN{N — l)/2. A finite hopping rate J lifts the degeneracy, 
the bands aquire finite widths and can even overlap if the tunneling parameter is large 
enough. 

The full energy spectrum calculated for N = L = 9 and J/U = 0.1, which 
corresponds to the Mott-insulator state, is shown in Fig. [5]^a). The lowest dot at K = 
is the ground state energy Eqi, see also Fig.[5^b). Also the lowest excited state is located 
at K = and has the energy Eq2- Together with the energies Ex2, where K ^ 0, they 
form the lowest excitation branch shown by the black lines in Figs. |5]^b,c). An increase 
of the system size leads to more dense distribution of the points and the solid line 
becomes smoother. We see that, at small momenta K, the lowest excitation branch can 
be approximated by a pseudo-relativistic form 

c^K = {Ek2 - Eoi)^ = (AS)^ + KW,s . (95) 

Thus we can estimate the effective velocity Vcs using the numerically calculated values 
of £^01) -^02 and Ell. The results for different system sizes are shown in Fig. [5t^d). With 
the increase of J/U the energy bands become broader and the gap in the excitation 
spectrum £'02 — -^oi becomes smaller. 

The energy eigenvalues in the lowest band in Fig. [5] correspond to particle-hole 
excitations of the form pi /ifc^l^Mott) with the total momentum kp + kh = K. When 



n,p 
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discussing the ground-state properties or the dynamics after a quench in the previous 
Sections, it was sufficient to consider translationally invariant states with K = 0, i.e., 
kp = —kh = k, where k corresponds to the relative momentum. In the discussion of 
the Sauter-Schwinger analogue, we considered a spatially constant potential gradient 
and absorbed it into the Fourier coefficients Tk(t) via a Peierls transformation, finally 
arriving at the evolution equations fl59|) and f l60|) for particle and hole operators with 
kp = —kh = k. However, for arbitrary potentials V^, this is no longer possible in this 
simple form. In order to satisfy the equations of motion ( TT91I2T1) for the correlation 
functions for an arbitrary potential V^, one should employ the generalised evolution 
equations 

idt-V,-^\p, = -^Y.T,,(h, + ^p^j . (97) 

Here the particle and hole-operators are fixed up to a /c/i^p-independent phase. For the 
limiting case V^ — )■ we find from Eqs. f l96|) and f l97|) the following eigenfrequencies 

1 



.h _ _ 

h 

1 



-{JT.^+uj,^)^ (9^ 



<= - ^l^iJTk, - ^k;) , (99) 

where ojk are the eigenfrequencies defined in ( l3Ti) . These excitations define the lowest 
band with the energies 

Ek^ = ujI - < , (100) 

where kh + kp = K. The spectrum is depicted in Fig. [6] and the effective velocity reads 

VeS = lVJi^U-J) (101) 

which has for J/U = 0.1 the numerical value Vcs/U = 0.269. Including the second order 
corrections in 1/Z^ we have a slightly lower value Ves/U = 0.264, compare Fig. [5](d). 

10.3. Probability distribution of the occupation numbers 

We calculate the ground state and then the probabilities p(n^) to have n^ atoms at a 
lattice site which satisfy the normalisation condition 

N 

As in the previous Section, we consider the system with N = L. From Eqs. ([3]) and (jlj) 
it follows that in the limit J = we have p{n^) = Sn,,,!, and in the opposite limit U = 
the probabilities are given by the binomial distribution 

jYj / 1 \ "^ / 1 ^ N-n^ 



pw-(F37a^UJ V-N^ ■ « 
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Figure 5. (a) Full energy spectrum and (b) its lowest part for A^ = L = 9. 
(c) Lowest part of the spectrum for iV = L = 13. (d) Effective velocity in a one- 
dimensional lattice with n ~ I atom per site and J/U ~ 0.1. 







Figure 6. Boundaries of the lowest energy excitation band in the continuum limit 
(solid lines) and energy excitations for L = 9 (points) from Eq. (|100p in a one- 
dimensional lattice with n = 1 and J/U — 0.1. 



The result for A^ = 14 at arbitrary J/U and at zero temperature is shown in Fig. [71 One 
can clearly see that the probability to have three particles or more at one lattice site is 
very small which is consistent with the approximations used in the 1/Z expansion. 

The particle-number distribution at finite temperature is shown in Fig. Mjo). 
Comparison with the zero-temperature result for the same system size [Fig. |8](a)] 
indicates that temperature has stronger infiuence at smaller values of J/U. 
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Figure 7. Probabilities to have n^ = (red), 1 (green), 2 (blue), 3 (magenta) atoms 
at a lattice site in a one-dimensional lattice of i = 14 sites with n = 1 atom per site 
at zero temperature. 





Figure 8. Probabilities to have n^ = (red), 1 (green), 2 (blue), 3 (magenta) atoms 
at a site in a lattice with n = 1 atom per site, (a) and (b): one-dimensional lattice 
oi L = 9 sites; (c) and (d): two-dimensional lattice of 3 x 3 sites; (a) and (c): zero 
temperature; (b) and (d): T/U = 0.3. 



10.4- Two-point correlation functions 

In this Subsection, we consider two-point correlation functions which have been 
discussed in Section O Due to the translational invariance, the correlation functions 
depend on the distance s = |x^ — Xjy| and have the property Fo{fi,iy) = Fo{s) = 
Fci{L — s) in view of the periodic boundary conditions. 

First we consider the parity correlation function F(_i)n(s). Its dependence on J/U 
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Figure 9. Parity correlation [(a) and (b)] and density-density correlation [(c) and 
(d)] in a one-dimensional lattice of L = 14 sites with n = 1 atom per site at zero 
temperature, (a) and (c): dependence on J/U for s = 1 (red), 2 (green), 3 (blue), 
4 (magenta); (b) and (d): Dependence on s for J/U = 0.1. Due to the periodic 
boundary conditions correlation functions increase starting from s = 7. 



as well as on the distance s is shown in Fig. |9l This result is in a very good quantitative 
agreement with the DMRG-calculations [20]. Note that our definition of J differs by a 
factor Z from that used in [20]. Fig. [9] shows the number correlation function -F^(s). 

Fig. [To] shows matrix elements of the one-body density matrix (6^6^+^) as well as 
its momentum distribution (!38|l . In a finite lattice the quasi- momentum takes discrete 
values which are integer multiples of 27r/L. These allowed values are marked in Fig. [TOT a) 
by dots. The momentum distribution calculated in the first order of 1/Z is also shown 
for comparison. We observe good quantitative agreement. 



10.5. Particle-number distribution and correlation functions in 2D 

The whole procedure of exact numerical diagonalisation described for one-dimensional 
systems can be generalised to higher dimensions. We did numerical simulations for 
two-dimensional square lattices of 3 x 3 lattice sites with periodic boundary conditions. 
Exact calculations for square lattices of the size 4x4 and larger were not possible due to 
the problem with computer memory. Due to the fact that the size of the two-dimensional 
system is very small one can expect strong finite-size effects. However, as we will see 
later numerical calculations give qualitatively correct predictions valid for large systems. 
The probability distribution of the occupation numbers is shown in Figs. [8](c) and 
(d). It is very similar to the one- dimensional case. In a lattice of 3 x 3 sites with periodic 
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Figure 10. (a): Momentum distribution (p8|) in a one-dimensional lattice with 
n — 1 atom per site at J/U = 0.1 and zero temperature. The result obtained by 
exact diagonalisation for a finite lattice of L = 14 is shown by dots. The solid line is 
a calculation in the first order of 1/Z for an infinite lattice, (b): Correlation function 
{w b^+s) calculated by exact diagonalisation for a finite lattice with the same values 
of parameters as in (a) . 



boundary conditions the distance s takes only three values s = 0, 1, a/2 which makes 
the study of the long-range behavior of the two-point correlation functions practically 
impossible. Nevertheless some useful information can be obtained in the Mott-insulator 
phase where the correlations decay exponentially. Fig. [TTT a) shows the dependence of 
the parity correlation function on J/U. As in the one- dimensional case, -F(_i)n(s) has 
a maximum at a finite value of J, which is, however, not in the Mott phase. The 
results in Fig. [TT] are in a good agreement with those obtained in [2U] by DMRG 
and MPS calculations for large systems. Correlations {n^nj) and (&t&i/) are shown 
in Figs. [TTT b). [TlT c). respectively. They become stronger for increasing values of J /U . 



10.6. Time evolution after quench 

If the complete set of the eigenvalues and eigenstates of the Hamiltonian is known, the 
time evolution of an arbitrary initial state | "1/^(0)) can be calculated exactly without 
numerical integration, provided that the Hamiltonian does not depend explicitly on 
time. The initial state can be decomposed into the eigenstates of the Hamiltonian as 

■Dk 

1^(0)) = J^X^c^oli^fi) , CKn = {i^mK^) ■ (104) 

K C=l 

If the parameters of the Hamiltonian do not depend on time, the evolution of the initial 
state is given by 

I^W) = EE^^r(t)|n^r) (105) 



with 



K r=i 



Vt. 



CKr{t) = ^ CKnCKnrexp {-iExnt) 



(106) 



n=i 
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Figure 11. Parity correlation (a), density-density correlation (b) and elements of the 
one-body density matrix (c) in a two-dimensional lattice of 3 x 3 sites with n = 1 atom 
per site at zero temperature, s ~ 1 (red), \/2 (green). 



The time dependence of the expectation value of any operator O can be calculated as 
(d)(t) = Y, J2^nKr\d\nK'r')c*Mt)cKT'it) . (107) 

Kr K'V 

We will be dealing with operators O which have the property 

(n;^r|C'|ni^/r') = {nKv\0\nKV')5KK' ■ 
Then the expectation value averaged over the evolution time is given by 

1 



:i08) 



{O) 



lim 

t— >oo t 



{0){t')dt' 



■Dk 



= J2Y1 $^(n/<r|0|n/^r') Yl l^i^^l' C*KnrCKnr • (109) 

RTF' f7=l 

We study time evolution of the initial state with exactly one atom at each lattice site 
which is the ground state with ii' = of the Bose-Hubbard Hamiltonian in the limit 
J = 0. Since the Hamiltonian after quench preserves the translational invariance, the 
time evolution involves only states with K = 0. In Figs. [12] and [I3] we present numerical 
results for the particle-number distribution p{n^) and correlation function {bj^by) in one 
and two dimensions. The purpose of this study is to address the question of (quasi) 
equilibration versus thermalisation in closed quantum systems. Time evolution of the 
probabilities p(n^) for one- and two-dimensional systems is shown in Figs. fT2ra) . fT3l(a) . 
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Figure 12. Quench and (quasi) equilibration in a one-dimensional lattice of L = 9 
sites with n ~ I atom per site, (a) Time evolution of the probabilities to have ti^ = 
(red), 2 (blue) atoms at a lattice site after quench J/U = — >■ 0.1. (b) Probabilities 
to have n^ = (red), 2 (blue) atoms at a lattice site for J/U = 0.1 as a function of 
temperature. Straight horizontal lines in both panels show the values of probabilities 
averaged over the infinite evolution time, (c) Elements of the one-body density matrix 
with s = 1 (red), 2 (green), after quench J/U = — >■ 0.1. (d) Elements of the one-body 
density matrix with s = 1 (red), 2 (green), for J/U = 0.1 as a function of temperature. 



respectively. On large time scales they oscillate around the averaged values shown by 
straight horizontal lines. For the chosen value of J/U = 0.1, the behavior p(0) is almost 
indistinguishable from that of p{2). Figs. [T2r b). [TST b) show the dependence of the 
probabilities on the temperature. Averaged values of the probabilities correspond to 
the effective temperature which is slightly less than 0.15 U. 

The time dependence of the correlation functions {bjj^bu) presented in 
Figs. [T2r c). [T3l(c) displays the same oscillating character. In the one-dimensional sys- 
tem, the effective temperature corresponding to the averaged values of {u'Jjy) can be 
defined but appears to be higher than that of the probabilities p{n^). In contrast, in the 
two-dimensional case, the correlation functions {b'^ubu) cannot be described by a ther- 
mal state, see Fig. [TST d). The absence of effective temperature in the two-dimensional 
system is consistent with the result obtained within the 1/Z expansion in Section [6l 



10.7. Tilt in one dimension 

In this Section, we calculate the probability to create a particle-hole excitation due to 
the time-dependent tilt. Initial state of the system |'?/'(0)) is the ground state of the 
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Figure 13. Quench and (quasi) equilibration in a two-dimensional lattice of 3 x 3 
sites with n = 1 atom per site, (a) Probabilities to have n^ = (red), 2 (blue) 
atoms at a lattice site after quench J/U = —J- 0.1. (b) Probabilities to have ni = 
(red), 2 (blue) atoms at a lattice site for J/U = 0.1 as a function of temperature, (c) 
Elements of the one-body density matrix with s = 1 (red), \/2 (green) after quench 
J/U = — >■ 0.1. (d) Elements of the one-body density matrix with s = 1 (red), y/2 
(green) for J/ C/ = 0.1 as a function of temperature. 



Hamiltonian ([T]) for finite value of J/U in the Mott-insulator phase. During the time 
evolution the system is described by the Hamiltonian ( 152|) with the on-site energies 

V, = Eof{t/T)x, , (110) 

where function f{s) has similar form as ( l72l) 



f{s) 



cosh'^ 



— cosh 



1 — cosh 



'1111 



_ 5 

with < s < 5, such that /(O) = /(5) = and /(5/2) = 1. 

In contrast to all the previous numerical calculation we do not impose anymore 
periodic boundary conditions. Instead of that we consider the case when the particles 
cannot tunnel between the lattice sites n = I and n = L and perform calculations using 
the basis of Fock states f l87|) . Numerical integration of the Schrodinger equation is done 
by means of the fourth-order Runge-Kutta method. The accuracy of the numerical 
integration is controlled by the conservation of norm of the state \ijj{t)). 

The results for the excitation probability per unit time 



|(^(0)|V'(5r))r 



;ii2) 
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Figure 14. a) Numerical results for the excitation probability per unit time in 
logarithmic scales for A^ = L = 12 and J/U = 0.1. b) Analytical results ([75]) for 
the excitation probability per unit time where the expressions for the effective mass 
(fTTj) and the effective velocity ([70|) have been used. 



where 5r is the total evolution time, are shown in Fig. IHT a). At short evolution times, 
the excitation probability Pgxc has a power-law dependence on Ae which corresponds to 
the perturbative regime of the pair production. One should keep in mind that finite-size 
effects start to play an important role if the evolution time exceeds L/ves, where Ves 
is an effective velocity for the propagation of excitations discussed in Section 110.21 For 
L = A^ = 12 and J/U = 0.1 this leads to the requirement tU < 7. At much longer times 
r, the dynamics of a finite-size system will be adiabatic and the excitation probability 
will tend to zero in contrast to the infinite system, where the excitation probability 
remains finite in the limit r — ;■ oo as determined by Eq. (1751) . 

In Fig. [TW b) we show the results of the same calculations obtained in Section [8] 
in the first order of the 1/Z-expansion where corrections due to time-derivatives of T^ 
have been neglected, see Eq. fl671) . The excitation probability (11121) and the particle- hole 
creation rate are related via 

p^^^ ^ l-(^Mott|p(oo)|^Mott) ^ l-il-2{p%)f ^ 2ivMM ,(113) 

where (p^^) = J2k l/^k| /N and /3k is the Bogoliubov coefficient defined in equation 
( !73l) . We observe a very good qualitative agreement with the results of exact numerical 
calculations, although the latter give somewhat smaller values of Pcxc- 
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11. Fermi-Hubbard Model 

Now, after having studied the bosonic case, let us investigate the Fermi-Hubbard model 
[H [121 EH] . We shall find many similarities to the Bose-Hubbard model - but also crucial 
differences. The Hamiltonian is given by 

H = -^Y. T,udlsCu,s + UJ2 ^l^i • (114) 

The nomenclature is the same as in the bosonic case ([1]) but with an additional spin 
label s which can assume two values s =t or s =4- In the following, we consider the 
case of half-filling {nl + hjj) = 1 where half the particles are in the s =t state and the 
other have s =J,. Note that the total particle numbers W = ^ n^ and A^^'' = ^ nj; 
for each spin species are conserved separately [H, W] = [H, N-^] = 0. The creation and 
annihilation operators satisfy the fermionic anti-commutation relations 

Cu,a, clbj = Sf,Jab , {C^,a, C^,b} = |4,a. c],,;,| = . (115) 

The fermionic nature of the particles has important consequences. For example, let 
us estimate the expectation value of the hopping Hamiltonian Hj. Introducing the 
"coarse-grained" operator 



y='S2Tf,uCu,s , (116) 



we may write the expectation value of the tunnelling energy Hj per lattice site for one 
spin species s as — J {&^^J^^^ j \fZ . This expectation value can be interpreted as a scalar 
product of the two vectors c^g \^) and c^ ,, |\E') and hence it is bounded by 

Inserting ||c^,, |\E') |p = (\E'| cjj,,c^,s |\I/) = (\E'|ri^,s |\I/), we get the expectation value of 
the number operator n^^^. In contrast to the bosonic case, this operator is bounded 
and thus we find \\c^,s |^) || < 1- Furthermore, the operator c^., in (11161) obeys the 
same ant i- commutation relations (11151) and thus we find ||c^s |^) || < 1 in complete 
analogy. Consequently, the absolute value of the tunnelling energy per lattice site is 
below 2J/vZ, i.e., decreases for large Z. 

The above result implies that the interaction term oc U always dominates (except 
in the trivial case f/ = 0) in the limit Z — t- oo under consideration. Hence, we are in 
the strongly interacting Mott regime and do not find anything analogous to the Mott- 
superfiuid transition as in the bosonic case. Note that often [971 [98] ^ different Z-scaling 
is considered, where the hopping term scales with J/\/Z instead of J/Z as in (lllip . 
Using this J/y/Z scaling, one can study the transition from the Mott state to a metallic 
state which is supposed to occur at a critical value of J where - roughly speaking - 
the hopping term starts to dominate over the interaction term. However, this transition 
is not as well understood as the Mott-superfluid transition in the bosonic case. With 
our J/Z-scaling in (lllip . we study a different corner of the phase space where we can 
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address question such as tunnelling in tilted lattices and equilibration vs thermalisation 
etc. 

11.1. Symmetries and Degeneracy 

In addition to the usual invariances already known from the bosonic case, the Fermi- 
Hubbard model has some more symmetries. For example, the particle-hole symmetry 
<^i,s ^ <^t^,s and thus fif^^s = cj^^^c^.s -H- n^,s = c^,sc),^^ = 1 - fif^^s is no longer an effective 
approximate symmetry, but becomes exact (for the case of half- filling considered here). 
Furthermore, there is an effective S't/(2)-symmetry corresponding to the spin 
degrees of freedom. To specify this, let us introduce the effective spin operators 

^^ = ^ E ^la <b c,,b = \ (nj - ni) , (118) 

ab 

and analogously S^ = J2ab^i,a<b^,^,b/'^ as well as S"^ = J2ab^la<bC^^,b/2 where a^^'' 
are the usual Pauli spin matrices. These operators satisfy the usual spin, i.e., SU{2), 
commutation relations and the Fermi-Hubbard Hamiltonian ( 1114^ is invariant under 
global SU{2) rotations generated by the total spin operators Stot = Xlu '^m- 

In the case of zero hopping J = 0, this global SU{2) invariance even becomes a 
local symmetry, i.e., we may perform a spin rotation at each site without changing the 
energy. As a result, the ground state (at half filling) is highly degenerate for J = in 
contrast to the Bose-Hubbard model (at integer filling). This degeneracy can be lifted 
by an additional staggered magnetic field (see Appendix [18]) and is related to the spin 
modes which become arbitrarily soft for small J. In this limit J <^ U, their dynamics 
can be described by an effective Hamiltonian, which is basically the Heisenberg model 
2 P 

with an effective anti-ferromagnetic coupling constant of order 1/Z'^. This effective 
Hamiltonian describes the Fermi-Hubbard Hamiltonian (I114p for half-filling in the low- 
energy sub-space where we have one particle per site, but with a variable spin S^. 

In order to avoid complications such as frustration for the anti-ferromagnetic 
Heisenberg model (I119p . we assume a bipartite lattice - i.e., we can divide the total 
lattice into two sub- lattices A and B such that, for each site in /i G ^, all the 
neighbouring sites v belong to B and vice versa. In this case, the ground state of 
the Heisenberg model (I119P approaches the Neel state for large Z 

he.i = ^^hinlnlfi, (120) 

which is just the state with exactly one particle per site, but in alternating spin states, 
i.e., s =1 for yU G ^ and s ="[ for u G B. Note that nj; is the projector on the \1) 

state hj^ = jl-'-) (l-''! while n^ projects on the \0) state etc. As usual, this state (I120p 
breaks the original symmetry group of the Hamiltonian (I114p containing particle-hole 
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symmetry, SU{2) invariance, and translational symmetry, down to a sub-group, which 
includes invariance under a combined spin-flip and particle-hole exchange etc. 

Let us stress that the Neel state f ll20p is only the lowest-order approximation of the 
real ground state of the Heisenberg model flllOp . there are quantum spin fluctuations of 
order 0{1/Z). These quantum spin fluctuations do not vanish in the limit J — > since 
J only appears in the overall pre-factor in front of the Heisenberg Hamiltonian (11191) 
while the internal structure remains the same. Only after adding a suitable staggered 
magnetic field (see Appendix [18]), the Neel state (11201) is the exact unique ground state 
(for J — !> 0). Either way, in analogy to the bosonic case, we can now use this fully 
factorising state (11201) as the starting point for our 1/Z-expansion. 

12. Charge Modes 

Starting with the Neel state (I120p as the zeroth order in 1/Z, let us now derive the 
first-order corrections. To this end, let us consider the Heisenberg equations of motion 

idtc^s = - -^^Tf^^c^s + Uc^sh^s (121) 

idtfi^s = - idtn^, = _ ^T^^ (at/^, - cj,/^,) , (123) 

where s denotes the spin label opposite to s, i.e., either (s, s) = {'\,i) or (s, s) = (|, t)- 
If we now insert these evolution equations into the correlation functions (cj^aCj^bJ^/^a^iys), 
{c^^a^yhri^afiyi), {claCiybfif^an^i), and {c^^a^uhn^ioXiub)^ wc find that they form a closed set of 
equations to first order in 1/Z, where we can neglect three-point correlations 

+ -^Tf,^{claCubnf,an^-b)o - -^Tf,^{claCi,hnf,aKb)o , (124) 

where the expectation values {n^a)^ and {n^b)o ^^ '^^11 as those in the last line are taken 
in the zeroth-order Neel state (I120p . In complete analogy, we obtain for the remaining 
three correlators 

idt{claCvbn^,-an^b) = + ^{n^a)o J^ T^^{c\^Cyb{n^a + n^a)Kb) 
- -^{nub)o X] ^'^^i^iaCKbnf.ain^b + n^b)) 

~ U {c^^C^bflfianub) 
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as well as 



and finally 



+ U {c^aCubnfManub) 

+ -^Tf,^{ciaCubn^,an^-b)o - -^T^u{claCf,bn^^an^-b)o , (126) 

^dticl^c^bn^aU^b) = +^{n^a)o Yl T^^{c\^c^b{n^-a + n^a)Kb) 
--^{nub)o Y TuK{claCKbnf,a{n^-b + n^b)) 

+ -^T^.A^laCubn^an^i)o - —Tf,^{claCi,.bn^,an^i)o- (127) 

We observe that the spin structure is conserved in these equations, i.e., the four 
correlators containing c'^Ci,^ decouple from those with c'^Ciyi etc. Thus we can treat the 
four sectors separately. Let us focus on the correlators containing c^ .c^i and introduce 
the following short-hand notation: If yU G ^ and u E B, we denote the correlations 

by {^ii^H'^^^t'^•4) = fjiv^''^ and (cj^^c^^^/.t'^^t) = flv^" ^ etc. Inserting the zeroth-order 
Neel state (1120p . we find four trivial equations which fully decouple 

^^tflf^ = - Uflf- , 
^dtflf' = + UfX- , 

idtfie'=Q. (128) 

Thus, if these correlations vanish initially, they remain zero (to first order in 1/2'). 
Setting these correlations (11281) to zero, we get four pairs of coupled equations 

jp, fOAOfl _ I Z \^ T f IsOs 
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^dtflf' = - ^ E ^-4''" - ^flf" ■ (132) 

Again, since these equations do not have any non- vanishing source terms (to first order in 
1/Z), they can be set to zero if we start in an initially uncorrelated state (see Appendix 
[T8l) . Note that they would acquire non-zero source terms if we go away from half- filling. 
The positive and negative eigenfrequencies of these modes behave as 



^^^U^VU^IJ^, (133, 

Thus we have soft modes which scale as u^ ~ J'^/U for small J and hard modes u^ ~ U. 
These modes are important for making contact to the t-J model [99] which describes 
the low-energy excitations of the Fermi-Hubbard Hamiltonian (I114p for small J away 
from half-filling. However, at half-filling, we can set them to zero. After doing this, we 
are left with four coupled equations, which do have non-vanishing source terms 

^^tf,C'^ = ^ E {^M./- '" - T.uf,f' } + Uflt'- - ^T,. , (135) 

^^tflf^ = ^J2 {T.J'f'' - T..flt- ] - Uflf- + ^T,. , (136) 

^dtfif- = ^J2 {T.Jt''' - T.uflf^} . (137) 

Due to the source terms JT^^/Z, these modes will develop correlations if we slowly (or 
suddenly) switch on the hopping rate J, even if there are no correlations initially. The 
eigenfrequencies of these modes behave as 



c^k= Vf/2 + 4J2T2. (138) 

A similar dispersion relation can be derived from a mean-field approach [96]. In contrast 
to the bosonic case, the origin of the Brillouin zone at k = does not have minimum 
but actually maximum excitation energy Wk- The minimum is not a point but a hyper- 
surface where T^ = (or, more generally, T^ assumes its minimum). After Fourier 
transformation of (I134p - (ll37p we find that the equations of motion conserve a bilinear 
quantity, that is 

dt [(/k"'" - 1) /k"'" + /k"'"/k"°"] = . (139) 

This relation holds, as in the bosonic case, also for time-dependent J{t). 
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13. Mott State 

13.1. Ground state correlations 

In complete analogy to the bosonic case, we now imagine switching J adiabatically 
from zero (where all the charge fluctuations vanish) to a finite value. Thus we find the 
following non-zero ground-state correlations 

flBlB ^ _ fOAO.A = ^ V - Tl - — ^ pi(xM-X.)-k Q4QN 

J fiu,gro\ind J fj.u,gT:onnd pr / j 9 \ Ld J ' V "/ 

k 

flBOA ^ fOAlB ^ ^ >p ::^ i(xM-X.)-k /...N 

J^i/jground J/^!^,ground at / j ' \ J 

k ^ 

Somewhat similar to the Bose-Hubbard model, the symmetric combination fll40p scales 
with J^ for small J while the other (11411) starts linearly in J. Other correlators such as 
{^IlI^h) ^^^ t)e obtained from these expressions. For example, if // and z/ are in A, we 
find, using ri^-|- -|- n^-|- = 1 and fi^^ + n^ = 1 

13.2. Quantum depletion 

In the zeroth-order Neel state (I120p . we have {n^-^h^^) = 0. Thus this quantity 
{n^-^fi^]) measures the deviation from this zeroth-order Neel state (I120p due to quantum 
charge fluctuations. In order to calculate {fi^-^h^])., we also need some of the other 
sectors discussed after (I127p . Obviously, the correlators containing cj^^c^rt- behave in the 
same way as those with &^<^c^]^ after interchanging the sub-lattices A and B. Thus a 
completely analogous system of differential equations exists for the correlations of the 
form (c|^^c^-|- 77,^1 ?7,,^|) = g]ji^^^ etc. If we insert (I123p in order to calculate idtiji^^n^:^), 
we find that these two sectors are enough for deriving (n^^n^i). Assuming ^ E B for 
simplicity, we find 






K^/i 



.-> hjLL ^ rill .1 hjLL .1 h^Li I \ / 



Setting the correlations with vanishing source terms to zero, we find 

k k 

Thus, in the ground state, the quantum depletion reads 

{finsfi^is) = {ri^su^s) = J^^2\^- — ] ■ (145) 

k 

As one would expect, this quantity scales with J^ for small J . 
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Figure 15. Time-dependence of the quantum depletion, the nearest-neighbour 
correlation function /^^°^, and the next-to- nearest- neighbour correlation function 
fuu^'^ in three dimensions after a quench within the Mott phase from J/U = to 
J/U ~ 0.5 in comparison to their ground state values. 



13.3. Quench 

Now we consider a quantum quench, i.e., a sudden switch from J = to some finite 
value of J, where we find the following non-vanishing correlations 



/. 






-I 



OaOa 
IJ,u,q\icnch 



1 ^ 2J'Tl ^ ~ ""'/^'-^^ e^(x.-x.).k ^ (14Q) 



J ^u, quench v^J^, quench/ t\j / > k'-^ 9 



,i(xp-Xy)-k 



k 



cu^ 



2 V^ 



sin(a;kt) 



3i(x^-x^)-k 



CJk 



(147) 



Again, these correlations equilibrate to a quasi-stationary value, which is, however, not 
thermal. For some of these correlations, this quasi-stationary value lies even below the 
ground-state correlation, see Fig. [151 The probability to have two or zero particles at a 
site reads 



\^pus"^^s) 



quench 



\fl ^s^ ^s ) 



quench 



^Y^^J'^i 



1 — cos(a;kt) 



OJt. 



(148) 
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This quantity also equilibrates to a quasi-stationary value of order 1/Z. In analogy 
to the bosonic case, this quasi-stationary value could be explained by a small effective 
temperature - but this small effective temperature then does not work for the other 
observables, e.g., the correlations. 

13.4- Spin modes 

So far, we have considered expectations values such as {cj^^Cubnfj,an^b)y where - apart 
from the number operators n^a and n^i - one particle is annihilated at site i> and one is 
created at site fi. These operator combinations correspond to a change of the occupation 
numbers and are thus called charge modes. However, as already indicated in Section [TT| 
there are also other modes which leave the total occupation number of all lattice sites 
unchanged. Examples are (cj^s'^iJ-sclsC^s) or {n^a'^ub) or combinations thereof. Many of 
these combinations can be expressed in terms of the effective spin operators in flllSp via 
(S^^Si). As one would expect from our study of the Bose-Hubbard model, the evolution 
of these spin modes vanishes to first order in 1/Z 

dt{S;Si) = 0{1/Z') , (149) 

consistent with the Heisenberg Hamiltonian flll9p . In analogy to the (n^n^)-correlator 
in the bosonic case, one has to go to second order O^l/Z"^) in order to calculate these 
quantities. Fortunately, the charge modes discussed above do not couple to these spin 
modes to first order in 1/Z and hence we can omit them to this level of accuracy. 

14. Tilted Fermi-Hubbard Lattice 

Motivated by the bosonic case, we now study particle-hole pair creation via tunnelling 
in a tilted lattice. Again, we assume a spatially constant but possibly time-dependent 
force on the particles which acts on both spin species in the same way. Accordingly, we 
modify our Hamiltonian via 

H = ~Y. T,uclsCu,s + UY, nlni + J^ V,{nl + nj^) , (150) 

/.il/,S fj, fl 

where V^(t) = E(t) -x^ denotes the additional potential. Performing the same procedure 
as before, we find modified equations of motion for the charge modes 

^dtltf-' = + ^ E (T^Jlf-' - T..fie^) , (151) 

if) {OaIb — _i_ _ \^ (t f IbIb _ T fOAOA^ 

+ (f/ + K - V,)fl-'- - ^T,, , (152) 

if) f IsOa - _ £ \^ (t f Isls _ T fOAOA^ 
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(153) 
(154) 



K^^,U 



In complete analogy to the bosonic case it is possible to factorise the differential 
equations for the correlation functions. We define effective particle and hole operators 
such that we have for the correlations functions without source terms 



{vIaKb) = flf' 



{Ka^^^b) = fl 



{K,bKb) = flf^ 

pOaOs 

ChlsKA) = flf^ 

{pIap^.b) = fie- 



(pIaP^a) = fie\ 
{pIbKb) = flf'' . 
ChlBP^^B) = f°e- , 

{KaP^a) = f'f\ 

{pIaKa) = flf\ 



JflUl 



and for the correlation functions with source terms 

(pIbKa) = flf\ {pIbP^,b) = flf- . 

This allows us to effectively factorise the equations for the correlation functions 



idth^A 
idtPij„A 



J 



j:^> 



K=/=fl 



MK \ h^^A + Pk,A ) + ( Y + ^M ) Pt^,B 



U 



K. 



+ V^ p^^A 



B + Pk,B 



-f... 



V,A 



u 

2" 



■idth^^B = ( -^ + ^M 1 Vb 



(155) 
(156) 
(157) 
(158) 
(159) 
(160) 

(161) 
(162) 

(163) 

(164) 

(165) 
(166) 



Due to the initial conditions, we can set the operators /z^^s and p^^A to zero. After 
Fourier and Peierls transformation fl58|) . we find the symmetric form 



idtPk,B = + -^Pk,B - JTi^{t)K,A ' 
idth]^A = - 'W^^'A - JT\^if)P\<i,B ■ 



(167) 
(168) 



where the Tk(t) are time-dependent. Now the line of reasoning is analogous to the 
Bose-Hubbard model. Initially, the operators evolve according to 

Ka = e+^^^/^ik , (169) 

]3k,B = e"^^*/25k , (170) 
with (y4j,^p) = and {Bl^Bp) = (5k,p. At the end of the evolution, we find 

/ik,A = fakik + /3kfik) e+*^*/2 ^^71) 
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Pkb = (/3k^ - «k^k) e-'""*/' . (172) 

In contrast to the bosonic case (where |akp — |/3kP = 1); we have |akp + |/3kP = 1- This 
difference reffects the fermionic nature of the quasi-particles and will have consequences 
for the case of resonant hopping (see below). The number of created particle- hole pairs 
then yields the depletion 

k 

Note that the equations (11671) and (11681) are analogous to the Dirac equation in 1+1 
dimensions, if we consider a small effective electric field E. In this case, particle-hole 
pair creation will occur predominantly in the vicinity of those points in k-space, where 
Tk vanishes, i.e., where the energy gap u^ in (I138P assumes it minimum. Inserting 
k = ko + 5k with T^p = 0, we may approximate Tk(t) via 

Tk(t) ^ [6k + A(t)] ■ [VkTk]ko . (174) 

Inserting this approximation into the equations (I167p and (11681) . we get 

^a^- JfVkTklko ■ [5k + A{t)]a^^ ■ ( ^^ j . (175) 

This is precisely the same form as a Dirac equation in 1+1 space-time dimensions if we 
identify the effective speed of light via 

Ceff = JlVkTjko , (176) 

and the effective mass according to 

rriescls = -. (177) 

Note, however, that Ceg depends on ko in general, i.e., the analogy only works if we single 
out a specific direction. In contrast to the bosonic case, we do not find a full analogy 
valid for all k-directions, since the dispersion relation is not isotropic near the minimum 
in the fermionic case. Nevertheless, we may use the analogy to the 1+1 dimensional 
Dirac equation in order to estimate the pair creation probability via 

where we have assumed a slowly varying field E. This result should be relevant for 
the investigations of the dielectric break-down in the Fermi-Hubbard model, see, e.g., 
[IOOlfTon[T02] . 

15. Resonant Tunnelling 

In the previous Section, we have studied the case of small potential gradients, i.e., small 
effective electric fields V"^(t) = E(t) ■ x^, for which we have obtained a quantitative 
analogy to the Sauter-Schwinger effect, which describes tunnelling from the negative 
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continuum (i.e., the Dirac sea) to the positive continuum. Now let us investigate the 
case of strong potential gradients. In this case, the lattice structure becomes important 
and resonance effects play a role. For simplicity, we assume a small hopping rate J <^U 
where we can solve the equations for the charge modes (I15mi54|) via time- dependent 
perturbation theory. In this case, Eq. f ll52p simplifies to 

{tdt -U-K + V,) fie^ = -^T,, + 0{J') , (179) 

as the other terms JiT^j.^fl^^^'^ — Tt^^f^^^^) are of higher order in J. We see that this 
equation becomes resonant if V^ — V^ = f/, i.e., if the energy gained by tunnelling from 
site /i to site v compensates the gap U + 0{J'^). In this resonance case, the correlation 
grows linearly with time /!]^^^ = iJtT^^/Z + 0{,P). Of course, Eq. (I153p has the same 
structure, but becomes resonant for the opposite case V^ — Vy = —U . Either way, the 
other two correlators 

^dtfl^'^ = ^Y1 {T,Jlf- - T^ufie^) , (180) 

and similarly fj^y^'^, grow quadratically /^^^'^ oc .Pt^. 

The same perturbative analysis can be done for the bosonic case, if we start from 
equations fll9ti2ip . Alternatively, we could employ the equations (12411260 in Fourier space 

[,5^ _ f/ + 3JTk(t)] fl^ = - V2JT^{t)Ul' + /f + 1) , 

[,a^ + U- 3.m{t)] fl' = + v^JTk(t)(/^^ + /f + 1) , 



^^Ji' = ^^Jl' = V2JT^{t)Ul' - h 



I') 



where we have inserted the time-dependent hopping matrix Tk(t) after the Peierls 
transformation fISS]) . Since T^ is periodic in k (the k-space is a periodic repetition of the 
Brillouin zone), the time-dependent hopping matrices Tk(t) are oscillating harmonicall}|§| 
for constant potential gradients. Thus the above set of equations corresponds to 
parametric resonance and can be analysed with Floquet theory. For simplicity, let 
us assume that the Tk(t) behave after the Peierls transformation as 

Tk(t) = I (e^^'V + e-'^°*Xk) • (181) 



In order to solve equations f l24ll26|) we make the Floquet ansatz 

oo 

fl' = E e^(^+")^°*/„^^ , (182) 

n=— oo 

oo ^ 

fl' = fl'=Y. e^(^+")^»V^^ - \ , (183) 

n=— oo 
oo 
f21 ^ ^ ^ri.^n)Eotf21 (^84) 

n=— oo 

§ For non-interacting particles, this behaviour is the basis for the well-known Bloch oscillations. 
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where v denotes the Floquet exponent and the f^^ are discrete Fourier coefficients of 
the correlation functions /^^. In order to satisfy equations fl24tl26p . the functions f^'^ 
have to fulffil the hnear system of equations 

(185) 



Mf = 


= 0, 










where we defined the infinite-dimensional matrix 
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(186) 



/ fl2 /■ll p21 fl2 fll f2l fl2 fll f2l \T 

[■■■1 J-ij J-ii J-i: Jo 1 Jo ' /o ' /i 1 Ji iJi ' •■■/ ■ 



:i87) 



(188) 



The set of equations (I185P has only nontrivial solutions if the determinant of the infinite- 
dimensional matrix vanishes, that is 

A(i^) = Det(M) = . (189) 

The above relation determines the Floquet exponent u up to multiples of 2tt and it can 
be shown that u satisfies the equality |lU3j 



sin2(7rz/) = sin^ A(0) . 



(190) 



If the hopping rate is much smaller than the potential gradient, that is J <^ Eo, we may 
expand A(0) in powers of J/Eo- Using the matrix-identity 

Det(M) = exp(tr{lnM}) 

we find up to forth order in J/Eo 



(191) 



Sm TTZ/ 



sm 



+ 



ttU 

J^IXkl^TTf/ 



Z^Eo{El - f/2) '°^ \-E~o 



Z^sin'[%)El{El-U^f{AEl-U^) 
SttU {AE^ - 5E2f/2 + U^) cos ^^""^ 



Eo 



Eo {-19E^ + 76E^U^ - 33U^) sin 



27rU 



(192) 



Two cases need to be distinguished. In the first case, the right hand side of (11921) 
is between zero and unity, the Floquet exponent is real and the correlation functions 
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are bounded. In the second case, the right hand side of ( 11921) is bigger than unity or 
smaller than zero, the Floquet exponent acquires an imaginary part and the correlation 
functions grow exponentially, /^^ ~ expC^utEo), corresponding to a Floquet resonance. 

We find the first resonance being located at U = Eq with a width of AU = 
2(53z/)max-Eo = A\/2J\x-k\/Z. The second resonance is located at f/ = 2Eo + 
16J2|Xk|V(3^o^^) and has the width AU = 2(Si/)max^o = l2^/2J^\xk\yiZ^Eo). In 
principle, there are resonances when Eq/U adopts higher integer values but they become 
smaller for increasing Eq/U . 

In contrast, the correlation functions in the Fermi-Hubbard model do not grow 
exponentially. This distinction between the bosonic and the fermionic case can already 
be deduced from the difference of the conserved quantities fl28|) and fll39p . While the 
relation ( l28|) allows in principle arbitrary large values of the correlation functions, we 
find from (11391) that f-^J^ ^ cannot exceed unity and is therefore bounded. 

16. Conclusions & Outlook 

For the Bose and the Fermi-Hubbard model, we studied the quantum correlations 
analytically by using the hierarchy of correlations obtained via a formal expansion into 
powers of 1/Z. Starting deep in the Mott regime J/U — f with exactly one particle per 
lattice site, we derived the correlations in the ground state for a finite value of J and 
after a quantum quench (i.e., suddenly switching on J). From these correlations, we can 
also infer the quantum depletion, i.e., the probability of having zero ("holon") or two 
( "doublon" ) particles at a given lattice site. It turns out that these observables approach 
a quasi-equilibrium state some time after the quench - but this quasi-equilibrium state is 
not thermal. Furthermore, we derived the particle-hole ( "doublon-holon" ) pair creation 
probability via tunnelling in tilted lattices and found remarkable analogies to the Sauter- 
Schwinger effect (i.e., electron-positron pair creation out of the quantum vacuum by an 
external electric field) in the case of weak tilts. For strong tilts, one obtains resonant 
tunnelling reminiscent of the Bloch oscillations for non-interacting particles. 

For the Bose-Hubbard model, we also studied a quench from the Mott state to the 
super-fluid regime and calculated the growth of phase coherence. Going up to second 
order in 1/Z, we derived the correlations of the number and parity operators, again 
in the ground state and after a quench. For the Fermi-Hubbard model, we found that 
the spin and charge modes decouple to first order in 1/Z. The dynamics of the charge 
modes (particle- hole excitations) already contributes to first order in 1/Z whereas the 
time-evolution of the spin modes requires going up to the second order in 1/Z, similar 
to the number and parity correlations in the bosonic case. 

Comparing our analytical results to numerical simulations for bosons on finite- 
size lattices in one and two spatial dimensions, we found qualitative agreement. Thus, 
although our analytical approach is formally based on the limit Z ^ 1, we expect that 
our results apply - at least qualitatively - to lattices in three {Z = 6), two {Z = 4), 
or even one {Z = 2) spatial dimension. There are only a few properties which strongly 
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depend on the dimensionality of the system, one example being the maximum of the 
parity correlator well within the Mott regime, which occurs in one spatial dimension 
only. In view of the tremendous experimental progress regarding ultra-cold atoms in 
optical lattices, for example, most of our predictions should be testable experimentally. 
In this paper, we used the zero-temperature Mott phase as our initial state - but 
the presented method can easily be applied to other initial states. For example, finite 
initial temperatures can be incorporated as well because our approach is based on density 
matrices. Even at zero temperature, it should be interesting to study other initial states. 
In the bosonic case, we could start with U = (instead of J = 0), i.e., in the super-fluid 
phase, where we may use p° = |tt)„(«| with the coherent state b^ \a) = a \a) , see, 
e.g., IH]. In this way, an order parameter (b^) 7^ is introduced at the expense of 
making the total particle number ill-defined [N,Pf^] 7^ 0. As another possibility, one 
could assume non-integer filling (n^) ^ N, where one has a non-vanishing super-fluid 
component even for arbitrarily small J. For example, taking (ra^) < 1, the initial state 
would be p° = \ip) A'ipl with \tp) = a\0) + /3 |1). In these cases, the time-dependence of 
p^ will be non-trivial in general. In the fermionic case, an analogous initial state would 
be p^n = l'^)u('^l with I'lp) = Oi lO^O^'') + (3 ll"'"!-''), which could describe a Bose-Einstein 
condensate of Cooper-like pairs, which may be stabilised by an attractive interaction 
f/ < 0, for example. 

Acknowledgements 

The authors acknowledge valuable discussions with M. Vojta, A. Rosch, W. Hofstetter, 
and others (e.g., several members of the SFB-TR12). This work was supported by the 
DFG (SFB-TR12). 

17. Appendix: Derivation of the hierarchy 

In this Appendix, we derive the hierarchical set of equations for the correlation functions. 
The quantum evolution of the one-site density matrix can be derived by tracing von 
Neumann's equation ([6]) over all lattice sites but p, and exploiting the invariance of the 
trace under cyclic permutations 

i^tp^, = ^ tr^ <^ ^ £„^p + ^ £f^p f + ^^<^ "^ 5Z ^"^ + ^^'P \ 

= \Y. ^l- trjp;. J + C,p, . (193) 

Using the definition of the two-point correlations given in ([8]), we arrive at ([9]). Similarly, 
the differential equation for the two-particle density matrix can be deduced by tracing 
over all lattice sites but p and z/, 

^^tp^,u = i {dtp^^" + Pf^dtpu + pudtPt^) 
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+ -^^luP^^v + >C^p^^ + C^p^^ . (194) 

With the definitions (|8]) and the time-evolution for the single-site density matrix fll93p . 
we find for the two-point correlation functions (fTOj) . The equations ([9]) and (fTOj) preserve 
the hierarchy in time if initially p^ = 0{Z^) and ff'^" = 0{1/Z) holds. In order to derive 
the full hierarchy, we define the generating functional 



J'(d) = T{{af,}) = In 



tr<^p(V)(l^ + d^) 



(195) 



where p is the density matrix of the full lattice and 



m,n 



are arbitrary operators acting on the Hilbert spaces associated to the lattice sites p with 
the local basis {|?7,) }. The role of this functional is to generate all correlated density 
matrices via the derivatives with respect to these operators d^ which are defined via 

If we consider an ensemble S = {pi, . . . , pi} of i different lattice sites pij^ . . . j^ pi, we 
obtain the correlation operators via 



d d d ^ 



(198) 

These operators are related to the corresponding reduced density matrix operator ps 
through the relation 

VJiPr=S i 

where the sum runs over all possible segmentations of the subset S into partitions Vi 
starting from the whole subset V = S and ranging to single lattice sites Vi = {p} 
where p'^^Lt ^ = p^ is understood. For two and three lattice sites, the above equation 
reproduces Eq. (|8]). 

Our derivation is based on the following scaling hierarchy of correlations: 

p%= O (Z^-l-^l) (200) 

where \S\ is the number d. of lattice sites in the set S. From the Liouville equation ([6]), 
the temporal evolution of J-" is given by 

idtJ'ia) = ^tr^ |d^£^^| (201) 

92 J^ dJ' d7 






tr„i. < («„ + au + aa&u) Cn 
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By taking successive derivatives and using the generalised Leibniz rule 

vuv=s 



d d 
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da 



as well as the the property 
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(202) 
(203) 

(204) 



we establish the following set of equations for the correlated density matrices: 



^dtPT 
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^J.es 



fJ.,u£S 






VUp=S\{fj.,u} 



'S <^corr "corr 
■'IikP{ij.}UV P{k}UV 



P{u}uv 



(205) 



^ ^ Z-^ E \^^"' Pft^y^'P PT^}yjv 

r QUQ=P 

— tr r^ /-corr _ , V^ -corr -corr _\ 

^'■u ^fiAP{f_L,u}UV ^ Z^ P{tJ.}UQP{u}UQ) 

'- QCP 

For i = 1 and £ = 2 we recover the equations (Q and (TTOj) . A careful inspection of 
this set of equations shows that the hierarchy in (I200p is preserved in time: Imposing 
the scaling p^™'' = 0{Z^~^^^) on the r.h.s. of the above equation, we find that the time 
derivative on the l.h.s. does also satisfy the hierarchy (I200p . Therefore, inserting (I200p 
into (I205p and taking the limit Z — )■ oo, we obtain the leading-order contributions 

Ate5 K^s ^J.es '- vcs\{ti} 

VUp=S\{fj.,u} 



idtpT 



E ^^pt 



r'S '-corr <^corr 

'~'^l^^P{^l}yJV P{k}\jP 



+ 



4e e 



QUQ=P 



n "Corr <-corr 

'-y^^ P{ii}yjv P{u}uv 



tr„ 



r'S \ ^ ;^corr ':corr 



}UC 



QCP 



PMup + (^(^^ 



-|5h 



(206) 



For i = 1 and £ = 2, we recover equations (I 111) and (TT^ . 

In contrast to the exact expression (I205p . the approximated leading-order equations 
(I206p form a closed set. The exact time evolution (I205P of the |iS|-point correlator dtffg" 
also depends on the higher-order correlation term ff^"^ involving |i5| + 1 points. The 
approximated expression (I206p . on the other hand, only contains correlators of the same 
or lower rank. This facilitates the iterative solution of the problem sketched in Section [3J 
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First one solves the zeroth-order equation fllip for p*^. Inserting this result p^ into the 
first-order (in 1/Z) equation flT2l) for p™'"'', we obtain a first-order result for p™''''. This 
first-order result for p^°J^ can then be inserted into the equation for p^™J which is of 
second order 1/Z^. Furthermore, we may use the first-order result for p™" in order to 
obtain a better approximation for the on-site density matrix p]^ which is valid to first 
order in 1/Z and contains the quantum depletion etc. Repeating this iteration, we may 
successively "climb up" to higher and higher orders in 1/Z. 

18. Appendix: Staggered Magnetic Field 

We assumed in Section [12] that the initial state of the Fermi-Hubbard system is given by 
the Neel state. However, for J = we have infinitely many states with same energy and 
vanishing total spin. In order to single out the Neel state, we add a staggered magnetic 
field to the Hubbard Hamiltonian, 

+ E (^^i^i - ^M^^i - K^^l) ■ (207) 

If we choose the magnetic field as A^^{x^ ^ A) = a, yl^|(x^ G i3) = 0, A^^{xy, E B) = a, 
and A^^{Xfj, E A) = 0, the Neel state is the unique ground state for J = at half filling. 
The Heisenberg equations f ll21l) - fll23p read now 

^^tc,, = - ^ E ^/^-^- + ^^A^^^M^- - ^^^^^^ (208) 

'Hs = ^ E ^MK^L - Ucln,, + A^,cl (209) 

idtfi^s = - idtUf^s = ^^Tf,^ (cLcms - cj,,c«,) , (210) 

To first order in 1/Z, we find the closed set of differential equations, cf. Eqs. (I124p - (ll27p . 

J 



idt{claCubn^an^b) = ^ ^ Tf,^{ciaCub{nf,a + nna)n^b){nfia)o 



-^ XI '^•^'^(^iaCKbnf,a{n^-b + n^b)){n^b)o 

z 



-^T'p^vLclj^vbri^afiyi)^ - —T^u{cl^c^bn^,anyi)Q (211) 

^dt{c\aCubn^,aKb) = 7 5Z T^^^^{^ia^vbifl^a + n^a)Kb){'n^,a)o 
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+ ^T^^{cl^c^bh^,an^i)Q - -7^T^u{claC^bnf,an^-b)o (212) 






+ {U + A^a - A^i,){c\^Cubn^,-aKb) 



+ -^T^y{cl^Cybn^an^-b)o - -^T^y{c\^c^bn^,-anyb)o (213) 






+ [Af^a — Ai,i,){c^^Cyb'n'naXI'yb) 
y-''fj,u\C^a^ub1T'fia1T'iyb/0 y 



+ -^T^y{cl^Cybn^any-b)o - -^T^y{c\^c^bn^,-an^-b)o ■ (214) 



If x^ G A and Xy E B we denote the correlations (cj^j^c^i '^/xt'^i/t) ~ fjiu^^^ 
(^M^^M-l-^Mt^i^t) — /°v^'^5 6tc. Inserting the zeroth-order Neel state, we find four equations 
which fully decouple, cf. Eq. (11281) . 

^dtflf^= -{U-a)flf- (215) 

^dtfX' = iU-a)flf- (216) 

^^tflf^ = (217) 

idtf,C'^=0. (218) 

In general, these four correlations are sources in the following four pairs of coupled 
equations, cf. Eqs. flT29|l -f lT32|l . 

^dtf,f^ = ^ E ^- ifj'^ + f.f') + <f^ (219) 

^dJif- = ^ 5^ T,. {f.f^ + flf^) - Uflf- , (220) 

^^tflf^ = - ^ E ^- Ulf"" + flf') - <f'' (221) 

^^*/i"'" = ^ E ^M. (/-'" + /- '') - <"'" (223) 
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^dtf,^'-' =^J2^>^^ (/- '" + /- '") + ^*'" ' (224) 

^dtfie^ = - ^ E ^- ifT" + fl^'") + <"'" (225) 

^dtflf' = - ^ E ^- if'^f" + Z^;^'") - ^/i'°" • (226) 

The eigenfrequencies of these equations read now 

, _ U + a±^U^Ti + {U-ay 
^k 2 ■ ^ ' 

In contrast to the eigenmodes (11331) we see that the modes fl227p have a gap in the hmit 
J — )■ 0. This enables us to switch on J and then switch off a adiabatically in order 
to have a well-defined initial state without correlations. Furthermore we have the four 
coupled equations, cf. Eqs. f ll34p -( 1T371) . 

^dtCf-^ = ^ E {^- (/°f ^ + fy^) - T.. iCf^ + f,f^) } (228) 

+ {U + a)fl-^--^T,, (229) 

-{U + a)flf- + ^T,, (230) 

^^tflf- = ^Y. {^M« (/- '" + /- '") - ^- (/m"°" + fl"'") } ■ (231) 
The eigenfrequencies of these equations read now 



w± = ±JU^Tl + {U + ay . (232) 



19. Appendix: Second-order Equations 

The differential equation for the three-point correlator reads 

„•/) "corr _ r- -corr , -*- r-S / - -corr , - -corr\ 
'''-'tPap-i — *-aPaP'y "r nr*-ap \PaP(3y "r PfSPay ) 



I J_ V^ ■ I pS ( fforr ^ I -corr -corr I 'corr -corr I ' ACorr\ 1 
+ r^ 2_^ ^J^K ^t-aK yPaMP'^ ^ P<^P Phi + Pa-t PkI3 + PaPnti-y) J 

' j.„ / /"S scorr" i /••^ A scorr \ 
—Pa^la ^^alsPpy Pa "T J-'apPfSPa-y j 
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~ -^Patr^ |L„^p^^ Pa + L^^p^p^p > 

+ (a -)■ /3, /3 -)■ 7, 7 -)■ a) + (a -)■ 7, 7 -)■ /3, /3 -)■ a) 

+ 0{l/Z^) . (233) 

In the following we use the matrix elements of ff°^^ and p^ in order 1/Z and define 

pTm= E P«"^'n«)«(«'l®l&)^(^'l®|c)7(c'l- (234) 

a,a' ,b,b\c,c' 

All three-point correlations can be deduced by permutation and complex conjugation 
from the following set of differential equations 

•Q ^001001 _ J sr^ rp / 001001 , ^7:^^00210 A 



J 
J 

'z 

J 



— \^ T fn'^owoi , /2„ooioi2\ 

E ^-/«7 (/S + v^/5) (235) 



J 



•Q ^001012 _ r r ^001012 I -^ V^ rp /^ ^001012 , ,/7j 002 

^<^tPap-y - -^Pap^ + ^ Z^ -'/3K(^PaK7 + ^ ^Pan^ 

■ ^ E ^- (^«r^ + 2pS") 

^ E ^««/S (/5 + V^/5) + ^r«7V^/5 (236) 



•Q ^002101 _ rr^002101 ^ V^ rp Z', /7T ^001001 , o^O 
«CtPa/37 -^Pafi-^ ~ ^ 1^ ^ 13k [^V ^PaK-y + ^Pc 



S^ 1^ /„002101 , ,/7J 002112 

^ Z^ ^7« (^Pa/3K + V2p„^^ 

K7^a,/3,7 

t^ X ^ ™ „io / „ii /r (.12 



!^a,/3,7 

E r.«/«7 (/5 + v^/3) - ^^"^v^/«7 (237) 



«;^a,/3,7 
J 



•Q ^002112 _ -' V^ rjn /^ AT 001012 , 9 00 

lOtPa/S-r ~ ~ ^ 1^ J^^Kyy^PaK-y + ^Pan 



002112 

K7^a,/3,7 
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- ^T,^ 72/21 ^ ^T„^y2/^2 (238) 



■o ^221001 _ "^ V^ T- / 221001 I . AT 22210A 
*C'tPa/37 ~ ^ Z^ ^ /^-^ l^P"«;7 + ^ ^Pa«7 j 

K7^a,/3,7 

\^ T Z' ^221001 I ^/;j^2210 
'y 2-^ TaKJaP [Jk-/ + V 2/^ 



J 



221012 \ 

re^a,/3,7 

K7 

K^a,/3,7 



r„,V2/S + ^T,/3V^/<i? (239) 



•o 221012 _ rr 221012 , ^ V^ rjn / 221012 , ^ A^ 2 

^OtPa/S-y - - f^ Pa/37 + 2^ 2^ ^/9k l^PaK7 + V 2p„ 



•^ X^ ™ / ooinn /:r 222112 

^aK7 
K7^o,/9,7 

^ $: T,. (v^p-r + 2psr ) 

I 5^ T..V2/S (/3 + V2/3) + ^T./3V2/^? (240) 



■Q 222101 _ rr^222101 '^ ST^ rp ( f^ 221001 , c) 2 

iOtPafS-r -^Po^M ~ ^ Z^ ^/3k l^v 2Pa«^ +2Pc 

K7^o,/9,7 

\^ T^ /^r,222101 I ,777 22211 
- ^ 2^ -'7«^ [PaPK + V2p^^^ 

«:7^o,/3,7 

^ E ^"'^v^/" (/S + v^/S) - ^^.7V^/S (241) 



Z 
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•n 222112 _ "^ X^ Ta /'A/9n221012 , n „. 



K7^a,/3,7 






■Q ^111001 _ "^ \^ T- /^ UllOOl , . /;^^ii2ioA 

^C'iPa/37 ~ ^ 2^ ^ /^-^ l^P"K7 + ^ 2Pqk7 ) 

K7^a,/3,7 

Y^ T^ /^lllOOl , ,/77^111012^ 



J 

'z 



f 

y - "7 J 0/3 y ""■ c«/^ ^ «7 

J 



if) nlll012 _ rr 111012 , ^ \^ T- /^n^llO^^ , ^/^„11^ 

^CiPa/37 ~ ~ ^ Pa/37 + ^^ ^ ^ /?« l^PaK7 +V^PaK 



12112 

K7 

«;7^Q,/3,7 



^ $: T,. (v/2pX°^ + 2pX") 



K7^a,/3,7 
J 

J 



'^a,/3,7 



•n ^112101 _ 7-r^ll2101 "^ \^ j, //l^ 111001 , 9„11 



412101 \ 

1K7 y 

«;7^Q,/3,7 
"^ \^ T^ /^/ill2101 I ^/;j^ll2112' 

«;^a,/3,7 



(242) 



-Ta-(faB 7^Tc,pf (243) 



K7^a,/3,7 
—^Tayfap -^Taisfa-y (244) 

J 
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J 



^—^07/0/3 "I — ^^0/3/07 (245) 



_L \^ T^ /^, /o ^112101 , 0^11211 

«;^a,/3,7 

- ^ E ^"'^ (/5 - v^/3) (/" + v^/" 

- ^T../^^ + ^T„,/- (246) 



■a ^200101 _ rr^200101 ST T /^n^OOlOl 4_ .fOr? 

^(^tPaPj - ^ Pa/37 ~ ~^ Z^ ^ /^"^ l,PaK7 + ^ ^Po 






201201 \ 

«K7 

K7^o,/3,7 
V^ T^ / ^200101 I ,777 200112 

;^ Z^ -'7« [PaUK + v2p^^«, 

K^«,/3,7 
J 



+ ^ Y. T.. {fj, - v^/^i) {f^l + v/2/,! 



K7^a,/3,7 



V2-^T. .11 V^-^T. .11 



•o ^200112 _ '^ \^ T f n'^OOlU , /^ 201212A 



j7 / J -^ pK yran-y ' * raK-y I 

Ky^a,l3,y 

, "^ \^ T^ A AJ ^200101 , 0^200112 
+ 2^ 2^ -'7« (^V 2p„;3« + 2p„^^ 

K7^a,/3,7 



K7^a,/3,7 
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+ -^-Ta-ffal3 -^Tafsfa^ (248) 

J 



•Q 201201 _ "^ V^ rp /^ 77^ 200101 , 9 201201 \ 



- — \^ T 1^^201201 , /^ 201212^ 
-^ 2_^ J 7K yPal3K + V ■^Pal3n J 

+ 4^T„,/- - 4^r.7/- (249) 



•Q 201212 _ rr 201212 , ^ \^ rp /^/7J 200112 , 9 20 

^"tPap-r - ~ ^Pq/37 + ^ Z^ J^I3k[VZP^^^ + ^Pan 



K7^Q,/3,7 



-L "^ V^ rp //7J 201201 I 9 201212 ^ 

+ 2^ 2^ -'7^ (^V2p^^«, +2p^^«, j 



+ ^ E ^"^ (/" - ^/") (/3 + v^/S 






■o ^310101 _ qrr „310101 ^ \^ rp / 310101 , ^A^^311201 
^CiPa/37 -"^^Pa^J ~ ^ Z^ ^/3k (^PaK7 +y^Pan-, 

k;7^q,/9,7 

\^ T^ / 310101 I ^/;j 310112 
-7 _^ -^ 7K I Pa/3K "T V ^Pq,/3k 

K^a,/3,7 
^ E •^S (•^«7 + ^/' 

"4^ E /"(/5 + v^/.^^ 

K7^a,/3,7 
J 



[•12 



^•-q ^310112 _ 9rr 310112 ^ Sr^ rp /^ 310112 , ,/^„3: 



311212 

K7 

K7^o,/9,7 
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^ 2_^ JaP ( /k7 + V 2/^^ 



Z 



E /"(/." + v^/3) - ^^«/^/" (252) 



.a,p^™ = 2t/p3-oi + ^ ^ T,. (v^P^ir + 2P^ 



N 7',3._ I A/'y'n-^""""" -I- ■/'/, 

aK7 

— V T /^n311201 , /^ 311212^ 

K^a,/J,7 

y 2_^ faP ( /k7 + ^'^Ik-i 

K7^a,/3,7 



z 



E /" (/3 + v^/3) - ^^^7/3 (253) 



•Q ^311212 _ rr 311212 , -' V^ rp //^ 310112 , 9 311212^ 

K^o,/3,7 

+ ^ E ^7. (v^pS""^ + 2pa'f ^) 

^^ E /"(/3 + v^/") (254) 



Z 

K7^q:,/3,7 



2 



By separating the two-point correlations in terms of order 1/Z and of order 1/Z 
according to p™" = p^™ + p'^fZ we find with 

pr^'^= E pr'""V)M(^'i ® i^).(^'i (255) 

m,m',n,n' 

the following set of differential equations 

■Q ^1001 ^ \^ T^ Z' ^1001 , ,77^2101 , ,A7 ^320 A 

^(^tp^,,y = ^ 2_^ ^t^K yp^^ + y^p^y + yip^u J 

\^ T^ f ^WOl I ,777^1012 I ,/q^1023^ 

- 2^ 2^ ^^« (^Pm« + ^2p^« + V3p^^ j 



~ "^ _^ -'a'kI Pooiooi + V 2p5q2ioi — PlllOOl 

— V/Pii2101 + V/P20OIOI + ^P201201 
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~ ~^ 2-^ -''^kI PlllOOl + V ^Piiioi2 ~ PoOlOOl 

— V^Po01012 ~ V /Po21010 ~ ^P021021 ) 

-3/o4^E(^^«/;^-^-/") (256) 

^^tpll'' = - ^ E ^- (^P^°" + 2P"" + v^Pf ") 

y 2_^ -^/^kI V ^^111012 T^ '^Pll2112 V -^^221012 '^P222112 

I AT ilKV |_ /Ti /^KI^ /^K!^ fry {IKV \ 

+ V0P310112 + VOP311212 " P2OOII2 ^ V^P201212 ) 
^ y _^ -'i^Kl V^P222101 + '^P222112 ^ ^ '^P\\2\m ~ ^Pll2112 

+ P022IIO + V ^Po22121 ~ V oPi32iio ~ VDP132121 ) 

- 3/0^ E ^^^^C - T.JII) (257) 

+ ^ E ^^'^ (^pi°°' + 2pi°." + v^pi^ 

~ y _^ -^ A«K I P0OIOI2 + V /Poo2112 ~ P1IIOI2 

— V^Pii2112 + V^P2ooil2 + ^P201212 
y 2.^ -^ i^K I V ^^^221001 ^ ^P221012 ^ ■^P\\\m\ ^Plll012 

I I//JK _, AT !^/lK AT I//1K AT V^1,K \ 

+ P02IOIO + V ^Po21021 ^ V oPi3ioio ^ vDPi3io21 ) 

- 3/0^ E (^^'^ (/m" + v^/") + v^^.'^ (/." + v^/")) 

+ 4/o^T^. (258) 
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^ J2 T,. {V2pT + 2plT + V6pZ'') 



J_ 



-f i/K I P0O2IOI + V /Po02112 ~ P112IOI 

— V^Pii2112 + v2pQ22iio + 2Po22121 ) 
+ ~^ 2^ -'mkI V2p221001 + ^P222101 " V ^PlllOOl ^ ^PlV. 



^311201 



T^ P2OOIOI T^ V ^P201201 V JPSIOIOI V U^3 






a/2J, 



-4/o^T^, (259) 



td,C^' = - 2Upll^' + ^ E ^- {pT + V2pfJ') 

~^ "^ 2_^ ^^1^ ( ^V OP22100I ~ V DP221012 + ^ ^PlZlOlO + ^P 



131021 



"... + fn- 



+ /o^T,. + /o4^ $: T..(/;i + v^/-) (260) 



■a ^3201 _ 9rr^3201 ^ V^ m /^3201 , ,/7J 3212A 
^CiP^i^ - 2tVp^^ - ^ Z^ ^^« (^Pm« + ^^Pmk j 



J 



J 



"^ _^ -'a'k I V OP2210OI + V DP22210I ~ V^Psioioi ~ '^P311201 



/3J 

Z '^'^ "" Z 



/o^T,. - /o^ 5: T,Ml + ^a (261) 



Kf^fJ.,^ 



~ "^ 2_^ p-i^ ( V 0P221012 + V DP222112 ~ ^ ^PziQii2 ~ 2p; 



311212 



'z~ 



/o^ E ^-(/.^^ + v^/.^D (262) 



K^H^U 



^^tpll'' = - UpfJ' - ^ E ^^^ {^pT + 2P"'^) 



Kf^fljh' 
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V 'JP222IOI ~ V OP222II2 + V ^Pl32110 + ^Pl32121 

(263) 



22\ 



Kf^fJ,,U 



^^,pll'' 



■J\^rp /^ ^101101 , ,777 211101 ,777^101112 9 211112 






- Pir°^ - V2pi°-- + v^p-;,i°i + 2p, 



^101112 

flUK 



ry / , ^ VK yP^KU 



111001 



_ 111001 

P^JLVK 



+ 



^2v^T,.(/;^ - /^ 



v^Pil.l^°^ - v^P^" - 2P. 
- V2piL^:" + ^^P'^^r + 2P, 



^112112 

fj,nu 



^112112 

flUK 



^<9tP 



2222 



Z 






(264) 



J 
J 



/^K 



^, /o ^102212 , 9^212212 ,/7J 512201 9^212212^ 
^V2p^^^ +2Pk^^ - V2p^^«, -2P^r.K j 



7/^ ^0000 



\^ T /^, /o ^221012 , 9 222112 ,/7J 222101 9^: 

- 2^ 2^ -'j-'c (^V ip^^,, + Ip^^y - V2p^^^ - 2p 



'lir) (265) 






--Tt 



Pk 



_ 001001 _ /^ 002101 I 001001 



|100012^ 



tP 



,0011 






K 

\^ 1^ / ^001001 , ,/7J 002101 ,/7T^001012 9^0021 
- 2^ 2^ -f^-K I^P^K^ + V 2p^«;, - Vip^^j, - 2p^^^ 



_ 001001 

PfJLUK 



-^^T,Mil-f, 



-v^p™+pr.r+2p; 



,002112 

flUK 



■21^ 



(267) 



z9jP 



1122 



- ^ E^-(/^^r' + v^p'r' - v^p^."^ - 2p; 



212212 

KUfJL 



_ 102201 

PfJLVK, 



^212212^ 



- V2p;r " + v^pjif °^ + ^P^ ^ 



- V2—T (f^^ - f 



fiuJ 



(268) 
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^dtPlT = - ^ E^-( - pTr - ^p" + p" + ^pIT') 



El^ ( /o„001012 I r) 002112 ,/7^ 002101 r, 002112A 
J'VK\y2,p^^^ + 2p^^^ -V2p^^^ - 2p^^^ j 



^dtfi.' = 


+ iu- 


-3JT^)f^'-V2JTM^' + fl' 


^dJl' = 


-iu- 


- 3JTk)/f + V2JUfi,' + fl' 


^dtfi.' = 


^dJl' = 


= V2JUf^' - fl') . 



+ ^'^T,^{C - fll) (269) 

19.1. Renormalised frequencies 

The two-point correlations to first order in 1/Z are determined by the differential 
equations 

source term, (270) 

source term, (271) 

(272) 

The 1/Z^-contribution of the correlations /k^,/k^,/k^ and f^ can be deduced from 
(I256II2591) . Defining the Fourier transform 

mm'nn' "}" \ "* mm'nn' ik-(x,j— x^) ('97Q'\ 

k 

we find from equations (I256H259P 

^^A''' = + upr - ^JUpr' - 3/o/^^) 

- V2.m {pr + pr - 3/o/^^ - 3/o/k^^) 

+ source terms , 

^^,p^ = - upr + ^JUpr - 3/0/f ) 

3/0/^^ - 3/0/f ) 

^^tp^ = ^d,pi''' = V2JT, {pi^^^ - pr - 3/0/^^ + 3/0/f ) 



+ source 


terms , 


- upr 


+ ^JUp'^'' 


+V2JT1 


c {pr + pr 


+ source terms , 


^d,pi''' = 


^ V2JT, {pi^^^ 


+ source 


terms . 



(274) 



(275) 



(276) 



As a next step we add equations flTTD]) and (1271) . (ITTT]) and (12731) . ( 1272]) and ( 1275]) and 
define 

pr = fl' + P^'' (277) 

pJc°^^ = /k'^ + P'r (278) 

pT^^ = fl' + P^^^^ (279) 

pr^ = fi' + pr^ • (280) 

From this follows a system of differential equations which is valid up to (9(1/Z^), 

^^tP^' =+[U- SJT^ (1 - 3/o)] pi''' - V2JT^ (1 - 3/o) {p^ + pI''') 
+ source terms , (281) 

^d,pi''' = -[U- 3 JTk(l - 3/o)] p'r + y2JTk(l - 3/o) {p^ + pl^ 
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+ source terms , (282) 

^^,p^ = ^^rp^ = V2JU1 - 3/0) {pr - pD 

+ source terms . (283) 

The homogeneous part of equations (I270p - (l272p is related to the homogeneous part of 
(I28ip -( l283|) via the substitution T^ — t- Tk(l — 3/o) from which follows immediately the 
renormalised frequency (!78|) . 

19.2. Parity-parity and particle-number correlations 

The parity-parity and the particle-number correlations are determined in 0{1/Z'^) by 
the differential equations (I264tl269p . Since the right-hand sides of (I264ti269p involve 
three-point correlations, we have solve the equations (l235p -( !M6l) . The calculations can 
be simplified by observing that it is possible to express the right-hand sides of (1264112691) 
by total time-derivatives using (l235l) . (l238ll . ( l239l) . ( l242|l . (l243|) . (l246|) and (12701) -(l272l). 
We find the exact expressions 



n---- — \^ z' 111001 I 112112\ / iq-x +ip-x +ik-x,^ , jq-x„+Jk-Xp+jp-x,A 



,1111 _ 

k,p,q 



4 E ifi'fi' + f^'fp) e^^-^'^^-^^-''^^ (284) 



iV2 

p,q 



r, — \^ 222112 / iq-x +ip-x +ik-x,^ , iq-x^+Jk-x^+jp-x„\ 

Pf^" j\r2 Z^ ^kpq {^^ "T e ) 



^2222 _ 

~ Ar2 

k,p,q 



p,q 

\ ^ ,001001 /^iq-x +Jp-x +Jk-x^ 1^ iq-x +Jk-Xu+Jp-x \ 

;p 2_^ Pkpq l^e ' " + e ) 

i«,p,q 



„0000 
^A"^ JY2 

k,p,q 



p,q 

1 

.0011 



— 11 ^ \ ^ 111001 ip-x^+ik-x„+iq-x^ 

^M'' M2 Z^ ^kpq ^ 



ATS 

k,p,q 



_ -'- \^ / ,001001 I ,002112\ 
^2 Z^ iPkpq -t-Pkpq ) 



^001001 _j_ ^002112^ iq-x^+jk-x^+ip-x„ 



k,p,q 



j-. E ifi'fi' + /"/") e^^''-^'^^-^^-'^^^ (287) 



p,q 

^1122 _ 1 V^ /^221001 I „222112\ ^ip-x„+Jk-x„+iq-x„ 

Pf^i^ ~ ~ Jp Z^ l^kpq + Pkpq j e ^ 
k,p,q 



+ ^ E pj,^p2q^^2giq-x,+ik.x,+*p.x, ^288) 



k,p,q 
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n022 _ 1 V^ ^221001^ip-x„+Jk-x„+iq-x„ , ^ V^ ^002112^iq-x^+ik-x^+ip-x^ 
k,p,q k,p,q 

-^^E/p'/^'^^^^^''^'^''''""^ (289) 

p,q 



aa'bb'cc' \ "^ aa'bb'cc' ik-Xc,+*p-Xo+iq-x ^900^ 



where we defined the Fourier transforms 

aa'bi 
Pafi-y j^2 

k,p,q 

After solving the differential equations for the three-point-correlations and inserting the 
solutions in f l284p - fl289p we find the parity-parity and particle- number correlations which 
are given in Section [9l 
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